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The primary intent of this work is to investigate 
computer aided compensator design using classical 
freguency response technigues in conjunction with the 
technigues of mpdern mathematical programming. A 
computer program to perform the automated design task 
is presented. This particular algorithm is based on 
the constrained optimization technique introduced by 
M. J. Bex . 

In particular, the desired epen lcop freguency 
response is specified for a number of discrete 
freguency points over the frequency range of interest. 
Then the ninimization routine is used to vary the 
compensator parameters in such a manner as to minimize 
a cost functional based on the difference between the 
actual and desired open loop frequency response of the 
compensated system. To illustrate the algorithm 
several detailed examples using the program are 
presented . 
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I 



INTBQDUCTION 



Concurrent with the proliferation of the digital 
computer as a viable, indeed often indispensable, 
engineering tool, the development of state-space methods 
resulted in the application of the digital computer in the 
area of modern control theory over the past decade and a 
half. Initially this application was primarily cf an 
analysis nature, but as the theory of optimal control came 
into wider publication and acceptance, the role of the 
digital computer began to take on aspects of the design 
function. This situation is graphically illustrated by the 
myriad of industrial proprietary and university developed 
computer programs for computer aided control system design 
based upcn state variable and optimal control theories. 
However, as pointed out by Mitchell [11] comparatively 
little work has been done to extend and apply the 
computational power, speed, and accuracy of modern digital 
computers tc the proven classical theory and methods of 
control system design. After this essential abandonment of 
the classical methods (particularly by the academic and 
theoretical factions of the control engineering community) 
in favor of the modern theory and techniques, recent 
publications indicate a realization of the need that both 
approaches to system design be used in a complimentary 
rather than mutually exclusive manner. In particular the 
works of Coffey, MacFarlane, Mitchell, Page, Rosenbrock, 
Stear, [5,8,11,14] and others show a renewed interest in the 
tried and proven classical theory cf design coupled with the 
new techniques and approaches developed by proponents of 
modern control theory, as applicable to the development cf 
usable computer aided control system design algorithms. 
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Ccntrcl system design based on frequency domain 
specif ications, for linear time invariant parameter systems 
using classical methods, has traditionally been done 
essentially using trial and error techniques. That is, once 
the form and location of the compensator has been decided 
upcn, the design process reduces to one of iterative 
analysis until the specifications are met or it becomes 
obvious tc the designer that the compensator form chosen 
will not satisfy the specifications. When systems cf any 
complexity are involved this method quite cften transforms 
into an extremely tedious and time consuming process. 

The availability of a good analysis program coupled with 
a digital computer operating in an interactive mode, to 
prcvide minimal turn around time, can be used to somewhat 
minimize the tedium of this type of design procedure. By 
varying the compensator parameters the designer may then 
observe the effects upon system performance without having 
to expend considerable effort in repeated calculations. 
This design by iterative analysis may be continued allowing 
the computer to perform the calculations until the designer 
is satisfied with the response obtained. With the recent 
refinements in mathematical programming techniques, it is 
possible to relegate some of the design decisions directly 
to the computer program, thus even further minimizing the 
time and efforts required of the design engineer. The 
ultimate coal, which is of course at present a long way down 
the road, wculd be to input only the plant dynamics and the 
desired specif icat ions and have the computer perform the 
design process completely, without any interaction with the 
design engineer. 

In this thesis an investigation of the automation cf the 
design prccedure is undertaken. In particular, giver the 
desired open loop response in the frequency domain a 
minimization routine is used in conjunction with an analysis 
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alqcritha tc systematically vary the compensator parameters 
tc approximate the desired response in accordance with a 
performance measure which indicates the deviation cf the 
actual open loop frequency response from the desired open 
lccp frequency response. 



15 



II 



THE AUTOMATED DESIGN PRO BLEM 



A. GENEBAL AUTOMATED DESIGN PHILOSOPHY 



The word design, as used in many engineering texts of 
various disciplines, runs the gamut of describing a process 
dealing essentially with scaling already existing 
configurations and components to one dealing with a highly 
creative activity allowing a maximum amount of freedom in 
satisfying a set of requirements or specif ications. On a 
cay to cay basis the majority of the design work 
accomplished by engineers lies somewhere between these two 
extremes, confining elements of relatively routine, straight 
forward calculations with those of creativity and 
imagination. 

The overall scope of the design process in general is 
indeed sc formidable that it often defies complete 
description in other than general terms, much less 
guantif icaticn in precise mathematical language. This 
latter is unfortunately a necessity in achieving the 
automation of this process within the limitations of the 
presently available hardware and software accessable by the 
engineer. Largely for this reason, when considering the 
automatic design problem it must by necessity be viewed, for 
the near future at least, in terms of a limited scope. That 
is, such items as creativity, experience, and intuition, 
which are undoubtedly elements of the total design process 
as viewed ir its entirety, cannot be achieved by a computer 
as the device exists at the present time. * 
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0 D€ cf tfce key elements of any design process is 
decision making. In the overall design process the engineer 
makes decisions based upon the information available to him 
at the tine cf the design undertaking. Quite often these 
decisions are based on a comparison of numerical values of 
various elements of the design configuration. It is in this 
specific area of the overall design process that the modern 
day computer can serve as a significant tool for the 
engineer involved in design, particularly when a large 
number cf ccapariscns (expressable in numerical terms) are 
invclved. In other words, the computer can be used in a 
combination analysis and elementary decision making role to 
aid the engineer in the overall design procedure. Ey means 
of analysis Eumerical values may be obtained, forming a 
basis frcm which comparisons can be drawn. The analysis 
portion is generally essential for exact design and forms 
the basis for any optimization. The speed available in the 
digital computer allows many calculations and comparisons, 
over a wide range of possible solutions, to be performed 
which would be prohibitive if peformed by hand. 

The employment of the digital computer as an aid in the 
aggregate design process is perhaps as much an art as the 
design process itself. While the use cf computer aided 
design pregrams can free the engineer from many burdensome 
and time consuming calculations, the individual engineer 
must still provide tie framework in which to interpret the 
results in the hard light of physical significance. The 
engineer must also be capable and willing to define the 
problem in a form with which the computer can work. More 
often than net, this latter requires a much sharper, clearer 
definition (and subsequently understanding) of the problem 
than if the same problem were being designed long hand. 
Unlike the analysis problem, which generally possesses a 
specific answer, the design problem may not have a unique 
solution which further complicates any computer aided design 
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algorithm. Another aspect that should not be overlooked is 
the fact that when using a computer aided design program the 
engineer often recei.ves, if nothing else, an indication of 
the soluticcs which are not feasible. This in itself is 
often net a trivial result especially in dealing with 
systems of any complexity. 

Another important factor which should net be overlooked 
is that once an answer is obtained using a computer aided 
design algorithm seme type of sensitivity study will more 
than likely be necessary. This is largely due to the fact 
that component values in physical hardware cannot generally 
be controlled to the precision with which a computer does 
its calculations. Therefore, while a particular algorithm 
for use in automated design may generate what is indeed a 
theoretically plausible solution the implementation of such 
a solution may not be possible if it is extremely sensitive 
to the parameter values involved. 



B. SOME ASEECTS CF COMPOTEBIZING THE CLASSICAL DESIGK CF 
CCMTBOL SISTEMS 



Classically the purpose of the design problem in 
feedback control systems has been to make a given plant 
satisfy a set of performance specifications imposed upon the 
system. resign of linear, time invariant compensators for 
the control of linear systems has formed a major portion of 
the design efforts of conventional ser vomechanisms. In the 
process of investigating the automation of this design 
procedure, certain key elements were found to be generally 
necessary in implementing a basic practical computer 
algorithm for automated design, regardless of whether it is 
based on time or frequency domain approaches. These 
elements consisted of: 



18 



1. A complete description of the plant 
cr process to be controlled. 

2. Specifications describing the 

response desired from the system. 

3. Some type of criterion function or 

f erformance measure relating the actual 
o the desired response. 

4. The form or structure of the 

compensator or control to be 

implemented. 

These items appear in varying degrees as a common 
denominator among present day automated design algorithms 
[1,5*11,13] regardless of the specific methods employed or 
systems being considered. That is to say, these elements 
appear in algorithms designed for time domain, transform 
domain, stochastic processes, and nonlinear systems. 
Generally items 1. and 4. are employed in some fcrm of 
analysis algorithm while items 2.. and 3. form the basis for 
a minimization (cr maximization) routine that provides a 
foundation for some type of decision making within the 
algorithm. This idea of using a performance measure upon 
which to base seme type of decision making routine is a 
direct carryover from developments .in the computerized 
design of optimal control systems which relies heavily upon 
the idea of minimizing some criterion function. Unlike the 
technigues cf optimal control however, the cost functional 
as used in this context is not necessarily a function cf the 
states and the applied control, tut rather a measure of the 
deviation of the actual response from some desired response 
which has been specified by the designer. This is perhaps a 
fine distinction but one that must be recognized and kept in 
mind. One of the many conseguences of this is that there is 
no longer a strong motivation to restrict the performance 
measure tc a guadratic form in order that the second method 
of Liapunov may be used to prove that a stable system will 
result. Bather, the form of the cost function in 
computerized design based upon classical control theory is 
chosen tc have properties which aid in the convergence of 
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the numerical minimization routine used in the algorithm or 
which has seme inherent physical significance in the design 
process. 



C. A FCE MO 1A1ICN Of THE AUTOMATIC DESIGN EBCBLEM IN THE 
EBECUENCI DOMAIN 

In this thesis the automatic design problem based on 
theory developed for classical design methods using 
freguency domain procedures is investigated. Several 
approaches based on frequency domain techniques have been 
presented by various people to automate the control system 
design process [5,11,13,14,20]. These algorithms range from 
methods designed primarily for use with multivariable 
systems to specialized interpretations of freguency domain 
data such as vector frequency domain techniques and inverse 
Nyquist methods. The suitability of these various methods 
have been demonstrated by their respective authors and 
others working in the area of computer aided control system 
design. fihile each represents a step forward in the line of 
progress toward automating the design process using 
frequency domain techniques the suitability of each 
algorithm is dictated to some extent by the particular type 
of design problem being solved. Many of these methods 
reguire gradient calculations, and those that do not rely 
primarily on simplex nonlinear programming algorithms. The 
former methods restrict the form of the cost function to be 
nirimized in the sense that it must be differentiable and 
problems inherent with numerical differentiating algorithms 
are obviously involved. The latter methods involving 
simplex algorithms have difficulty in handling constrained 
minimization problems, particularly when solutions exist 
near or on the constraint boundaries. 
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One of the advantages of the frequency response approach 
lies in the fact that it lends itself easily to the use of 
open-loop frequency response data in designing closed loop 
systems. The use of frequency response methods provides a 
graphical portrayal of the system response that lends itself 
well to quick interpretation by the designer. Also since 
approaches such as those developed by Bode were intended to 
give simple approximate methods for use in hand 
calculations, adaptation of these methods in a computer 
algorithm allcws the engineer tc plan his design using such 
techniques and then employ the aid of the computer without 
the requirement cf significant problem restructuring in 
order to cast the problem in some particular format not 
generally familiar to the designer, but necessary for 
solution by a computer algorithm. With systems involving 
complex transfer functions such restructuring can itself 
become quite burdensome thus defeating one cf the primary 
reasons for turninq to computer aided design. Perhaps it 
should be pointed cut at this point that even using a 
computer algorithm to aid in the design process does not 
totally relieve the engineer of all hand calculations. Time 
must be taken to carefully map out a plan of attack on a 
particular problem. In essence the designer must throughly 
understand the problem and this is generally achieved by 
seme preliminary rough sketches and calculations before the 
computer ever makes its appearance in the design process. 

Throughout this thesis, the vehicle used for presenting 
the automated design algorithm is the single-input 
single-output unity feedback system of the form shown in 
figure II-1. 
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figure II- 1 

Ibe plant is considered to be linear, time invariant, and 
ccupletelj known. Thus the plant may te described by a 
rational transfer function in the Laplace variable s of the 
form 

n 



G(s) 



I I (s + z ) 
i=0 i 

fr (s + V 

i=0 



where z and p may he complex constants. The compensator 
i i 

G (s) is also considered to be a rational transfer function 
c 

in s. 



The specif icat ions are given in terms of the system open 
loop frequency response, and the desired response is 
achieved by reshaping the system open loop frequency 
response through the use of the compensator in the feed 
forward path. This basic procedure is well known to most 
control engineers and is by no means a novel approach to the 
design problem [ 3, .10, 17, 1 8 ], but through the years it has 
found extensive use in the practical design of control 
systems. 
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The method proposed and under investigation here is a 

relatively straight forward adaptation of this technique 

with a computer algorithm harnessed to do the compensator 

parameter variation in order to systematically reshape the 

system open loop frequency response with a minimum amount of 

calculation and iteration required on the part cf the 

designer. The desired open loop response is specified by 

means of a magnitude and phase profile ever the freguency 

range of interest. This is nothing more than a set of 

desired magnitude and corresponding phase values at discrete 

frequency points selected by the designer. In addition a 

specific form or structure is postulated for the compensator 

transfer function G (s) . At this point cf the automated 
c 

design process the algorithm relies heavily on the 

experience and intuition of the individual designer in 
choosing a structure or form of compensator that will have 
seme reasonable probability of success in reshaping the 
frequency response curves in order to satisfy the design 
specification. Since no firm analytical basis exists for 
the determination of a specific structure for the 

compensator, the computer algorithm must be given this 
information or an iterative method based on some type of 
criterion must be used in order to establish the form. This 
latter method was felt to be much too time consuming and 
restrictive, because of the wide range of variations in 
possible structures that might have to be considered, and 
essentially beyond the scope of this investigation. Thus 
the choice of the compensator form is relegated to the 
program user as a design decision that must be undertaken 
prior to using the algorithm to determine specific values 
for the coefficients of the compensator transfer function 
that will yield the desired response. This procedure should 
net present a majpr stumbling block for the designer 
provided that he possesss a thorough understanding cf the 
problem to be solved. In fact, by inputing the specific form 
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of the compensator the designer may investigate compensator 
configurations which might otherwise be considered 
unfeasible for the solution of the problem if some type of 
iterative algorithm were used in determining the structure. 
With these items forming the basic inputs to the algorithm a 
nonlinear performance measure coupled with a minimization 
routine is then used to minimize the difference between the 
desired and actual frequency response of the open loop 
transfer function gi.ven by, 

T (s) = G (s)»G (s) 
c 

The performance measure (or cost function) selected is of a 
normalized form primarily as a matter of convenience in 
interpreting the algorithm’s results regardless of the 
particular problem being considered. That is , the cost 
function has a general form given by 



A ( ja>) 

J (ja>) = f (1.0 ) 

D (3if) 

where D ( jo?) represents the desired open loop frequency 
response of the entire system, A(ja>) is the actual open loop 
response, ard w represents the set of discrete frequency 
values at which the cost function is to be evaluated. As 
mentioned previously this general form of normalized cost 
function mas selected primarily as a matter of convenience 
in that a particular cost function can be constructed in 
such a manner that the optimal or smallest value of J ( jw) 
under ideal scluticn of the problem will be zero. 

By inputing the desired response in the form of the 
magnitude and phase profile described above, normal 
frequency domain specifications such as gain margin, phase 
margin, and bandwidth (here defined as the open loop gain 
crossover frequency) may be supplied to the algorithm in one 
type of format. With these specifications forming part of 
the gait and phase profiles they are automatically 
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incorporated into the overall cost function each tine the 
algorithm evaluates the frequency response of the open loop 
system ever the range of frequencies specified ty the 
designer, thus also eliminating the requirement of 
specialized computation within the algorithm to check for 
satisfaction of these particular specifications as the 
compensator parameters are varied. 
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Ill 



THE M INIM IZATION TE CHNI QUE USED 



A. GEN EE A I 



A variety of functional minimiza tion techniques have 
come into use in recent years are discussed in ref. [ 8 ] - 
The minimization algorithm used in this work was developed 
originally ky M • J. Eox, [4], and is capable of finding the 
minimum cf a general nonlinear cost function, composed of 
several variables, within a constrained region. The complex 
method of Eox has teen implemented in program form as part 
of the standard subroutine library at the Naval Postgraduate 
School S. E. Church Computer Center. The subroutine 
(BCXPLX) , was originally programmed by fi. E- Hilleary cf the 
staff of the Naval Postgraduate School Computer Center and 
has subsequently been modified slightly by the author in 
order to render it more suitable for use in the scheme of 
the compensator optimization program. 

In addition to its availability as a standard 
subroutine, several other features of this particular 
minimization algorithm made it an attractive choice for use 
in minimizing a cost function based on frequency domain 
specif ications . iihile admittedly more efficient and 
sophisticated algorithms for function minimization are 
available, the flexibility and ease of programming of the 
complex method of Bex are among considerations that should 
not be overlooked. This method does net require that the 
derivatives cf the cost function or implicit constraints be 
calculated, thus avoiding the sundry problems typically 
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associated with numerical methods for the calculation of 
derivatives, and eliminating the restriction that the 
derivative of the cost function exist over the range of 
interest. In implementing the minimization method the cost 
function need not be expressed explicitly in terms of the 
variables tc be manipulated by the program in the 
minimization process. Unlike many of the other minimization 
programs available, BOXPLX, also provides a restart 
capability. That is, once a possible minimum value cf the 
cost functional is reached the program automatically 
restarts from racdpmly generated initial values cf the 
variables. fchile this method by no means guarantees ttat a 
global minimum will be achieved for every conceivable 
problem the restart capability does at least somewhat 
attempt tc avoid the problem cf the minimization terminating 
on a local minima within the search area. This method is 
also tc seme extent scale independent in that the size of 
the initial cecmetric figure used in the search for the 
minimum value cf the objective function is rcughly scaled to 
the order cf the problem variables, through the use of the 
difference between the lower and upper bounds on the 
variables. 

As is the case with other minimization algorithms, this 
particular ere is not void of inherent disadvantages and 
difficulties. Present in any minimization algorithm of 
course is the recurring problem that no guarantee can be 
made of achieving a global minimum for all classes of 
problems. The complex method cannot handle equality 
constraints without modification cf the algorithm. Also as 
the number of variables increases the method rapidly becomes 
inefficient. Finally the unconstrained problem is 
apparently mere efficiently handled via gradient techniques. 



B. DESCBIF TICN OF THE ALGORITHM 
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The n inimization method developed by Ecx has its basis 
in the sequential simplex techniques of linear programing, 
but is specifically designed to obtain solutions to problems 
cf constrained optimization and to avoid such probleas as 
the inability of constructing the first siaplex. Thus the 
resulting scJuticn of the minimization problem is a set of 
variables which nay lie at the extreme edges or boundaries 
cf their permissible ranges. The problems generally 
associated with constraints on the variables is essentially 
attacked by using a flexible geometric figure having n+1, or 
more, vertices and capable of expanding or contracting, in 
any or all directions. Therefore, unlike many of the 
siaplex algorithms no provision is made to maintain a 
regular cecaetry between the various points in an 
n-dimensicnal space, but rather a changing structure (dubbed 
a complex), capable of flattening, rounding corners, and 
eventually collapsing on the point representing the minimum 
of the objective function, is employed rendering the method 
more flexible in handling constrained minimization problems. 
The number cf independent variables in the objective 
function and constraint equations dictate either directly or 
implicitly tbe dimensionality cf the space to be searched. 
The vertices making up the complex are continually rejected 
and generated as a search for the minimum cf the objective 
function is carried out with each new vertex generated being 
reguired to satisfy all the imposed constraints. 

The complex method searches for the minimum of an 
objective function J (x) , where the vector x represents a set 
of n variables, within a region of space bounded by upper 
and lower Units or the variables (explicit constraints) 
given by, 

L < x < U , i=1 ,2,3,. . . ,n 
i i i 
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and in addition subject to tbe restrictions imposed b y the 
constraint functions (implicit constraints) cf the fora, 
h^(x) > 0 j=1,2,3,...,m. 

As with the simplex methods this is an iterative algorithm. 

It uses a set of k > n+1 points which form the vertices of 

the complex and simultaneously satisfy all the imposed 

constraints. Initially only one point or vertex must 

satisfy all the constraints and serves as a means of 

generating the other k-1 vertices through the use of 

pseudorandom numbers r , uniformly distributed over the 

i 

closed interval frpm zero to one. The additional k-1 

vertices cf the first complex are established by the 
following formula; 

x = L + r • (U - L ) 
i i i i i 

As can he seen from the formula, the vertices generated in 
this fashion will always satisfy the explicit constraints, 
however, each one must be checked to insure compliance with 
the implicit constraints,. If an implicit constraint 
violation is indeed found to occur at this stage then the 
trial vertex is simply moved, repeatedly if necessary, 
halfway in toward the centroid of the ether already accepted 
vertices, until ultimately a feasible point is found. 
Continued application of this procedure will result in the 
generation cf the k-1 additional points required for the 
initial complex within the feasible region defined bj the 
constraints. 

Cnee the initial complex has been determined the 

iteration proceeds by evaluating the objective function J (x) 
at each vertex. The vertex at which the objective function 
has the largest value is then designated as the current 
worst vertex and this point is reflected through the 
centroid of the remaining vertices establishing a new 
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complex. If this worst vertex is designated as x thee the 

overref lection is accomplished according to the formula 

N CR w 

x = (1.0 + alpha) «x - alpha«x 

CR 

where alpha > 1.0 is the overreflection coefficient, x 

N 

represents the centroid of the remaining vertices, and x is 

the new vertex. If the overreflection of the current werst 
vertex results in a violation cf the constraints or if the 
new point still has the worst value in the set of vertices 
then the pcint is moved halfway toward the centroid used in 
the over ref lection process. This retraction is repeated (if 
necessary) until the overreflection coefficient is reduced 
to some snail value, delta, or a suitable new vertex is 
established. If the overreflection coefficient is reduced 
to the value delta and the objective function value at this 
new point is still the worst in the set of vertices then the 
projected vertex is replaced by its original value and the 
second worst vertex will be overreflected instead. This 
process keeps the complex moving toward the minimum of the 
objective function provided the figure has -not collapsed 
into its centroid. Once assured that the new complex 
satisfies all the constraints and results in an improvement 
of the objective function, vertices of the complex are 
rejected and generated in a systematic manner aimed at 
minimization cf the objective function. Thus essentially 
the complex moves over the feasible region eventually 
straddling and collapsing upon the point that renders a 
mirimum value for the objective function. 

This process continues and is eventually terminated when 
the complex shrinks to a predetermined acceptable small 
value epsilcr given by: 
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< epsilon 



/ 1 r^\ C V 

{ k 2^ c J ( - } ’ J ( - } ] 

1=1 

c 

Where J (i ) represents the value of the ob jective function 

V 

at the centrcid and J (x ) is the value at the various 
vertices. 

As pointed out by Box [4], the exact values cf the 
overref lecticn coefficient, alpha, and the exact number of 
vertices, k, do not appear to be critical to the functioning 
of the algorithm provided that they are greater than unity 
and n+1 respectively. The main concern here is that a value 
of alpha larger than unity keeps the complex from shrinking 
prematurely as it traverses the feasible region in search of 
the minimum objective function value. Values of alpha = 1.3 
and k = 2n as given in ref. 4 are used here since the exact 
values choser are apparently not critical and good results 
have been reported, by the originator of the algorithm, 
using these particular values. For a detailed example of 
the iteration process and a more sophisticated explanation 
of the algorithm the interested reader is directed to refs. 
2, 4, and 6. 
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PEOGEAM DESCR IPT ION AND UTILIZATION 



IU. 



A. DESCRIPTION Of ALGORITHM IMPLEMENTATION 



The frequency domain algorithm described in the latter 
part of Chapter II has been implemented in program form 
using FORTRAN IV coding and configured for use on the IBM 
360-67 computer system available at the Naval Postgraduate 
School Computer Center. With the exception of the graphing 
subroutine the remainder of the program is self contained 
and may be used on any computer system having a standard 
FORTRAN IV compiler. Without the overhead routines 
reguired of the graphing subroutine, the basic program 
requires approximately 180k bytes of core memory, and 
including the core requirements of the plotting routines 
roughly 210k bytes of memory are required by the pregram. 
While admitedly less core would be necessary if such 
techniques as array overlaying were employed when 
programming the algorithm, it was felt that for simplicity 
in following the flow of the program by future users these 
techniques were best left unused. 

The computer aided linear time invariant compensator 
optimization program (CALICO) is essentially composed of 
four main parts with numerous additional small subroutines 
functioning in supportive roles. Since often a program user 
may wish tc perform certain modifications to programs of 
this nature in order that they may better serve his 
particular purpose. Appendix A provides a description cf the 
supportive subroutines used within the program to facilitate 
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accomplishment of any program mcdif icaticns that may be 
desired. A block diagram of the major components cf the 
overall program are shown in figure IV-1. The main program 
serves primarily tc input data describing the system under 
consider aticc and to output the values of the roots and the 
gain of the compensator that will best achieve the desired 
response in accordance with the cost function selected. 
During program execution, the main program also serves to 
direct the seguence of execution of the other major sections 
of the C A1ICC algorithm. 

Since the plant and compensator are assumed to be linear 
time-invariant, they may be completely described by transfer 
functions in the independent Laplace variable s. Once the 
plant and assumed fprm of the series compensation are input 
along with the profiles forming the desired phase and 
magnitude response program control proceeds tc the 
minimization subroutine, BOXPLX. From the point of view of 
computer aided control system design, associated with this 
subroutine are two key function sub-programs, FE(X) and 
KE (X) . The function sub-program FE (X) performs the 
calculation cf the system open loop frequency response, at 
the discrete frequencies provided in the input data, and 
calculates the value of the cost function based upon the 
difference between the actual and desired frequency response 
at the discrete frequency values. Six separate cost 
functions, from which the program user may select one, are 
incorporated into this function sub-program. The value of 
the cost function is used by the minimization routine to 
vary the compensator parameters in order to minimize this 
guantity. Ibe function sub-program KE (X) checks fer any 
violation c i the implicit constraints, which in this 
particular case are concerned mainly with the presence of 
right half plane roots. The implicit constraints used, also 
impose limits on the real and imaginary parts cf the 
numerator anc denominator polynomials evaluated at any 
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MAIN PROGRAM 



Read input data. 



Output optimized results; call simulation, 
stability, and plotting routines. 



SUBROUTINE BOXPLX 

Systematically varies the compensator 
parameters to minimize the cost function J. If 
the number of trials is exceedid the best 
solution achieved to that point is selected as 
the solution. If the cost function remains 
unchanged for a specific number of trials the 
minimum is considered as having been achieved. 



A 



y_ 

FUNCTION FE 
Calculates the fre- 
quency response of 
che open loop system 
at the discrete val- 
ues specified in the 
input and forms the 
difference between 
desired and actual 
response making up 
the cost function. 



T 



i__j 

FUNCTION KE 
Checks for right 
half plane roots 
and violation of 
other implicit 
constraints . 



SUBROUTINE SIMULA; Calculates the open loop 
frequency response using optimized parameters. 



- — ... - .. a — . — 

SUBROUTINE GRAPH; Plots the desired graphs 
specified bv the user. 



Program Block Diagram 
Figure IV-1 
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particular freguency. These latter limits are a direct 
result of the finite computer word size and the particular 
method csed to represent numerical values within the 
machine's memory. Beth of these factors have an effect in 
determining the maximum and minimum size numerical value 
that may he represented within any particular computer. 
While the ccefficient values determined by the minimization 
routine will net exceed these maximum or minimum values, the 
real and imaginary parts of the numerator and denominator 
polynomials evaluated at a particular freguency are sguared 
during the process of determining the magnitude cf the 
freguency response. Without the implicit constraints these 
sguared values can exceed the size of the numerical 
capability of the machine, particularly in the case cf high 
order systems, causing premature termination of the 
algorithm anc hence failure in obtaining a solution. Cnee 
the minimization subroutine determines the coefficients that 
result in a minimum value of the cost function or the 
minimization process is terminated for other reasons, 
program ccntrcl is transfered back to the main pregram. 
Here the roots and gain of the compensator which resulted in 
a minimum cost function are calculated and output. The open 
loop freguency response of the sytem is then calculated 
using a mere densely distributed number of freguency values 
than would normally be selected by the designer for use in 
the minimization process. This is done largely to show the 
existence cf any resonant peaks which might otherwise be 
straddled and not apparent if only the discrete freguencies 
selected ty the designer were used in the plotting routines. 
The stability cf the entire system is then checked by adding 
the open lccp numerator and denominator polynomials and 
applying tbe Rcuth stability criterion to the resulting 
polynomial which is the system characteristic eguation. 
Finally the results at the discrete freguencies specified by 
the designer and the "continuous” curve generated ty the 
SIFULA subroutine are output on Bode and gain versus phase 
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plots. In addition, plots of the difference between the 
desired and actual results at the discrete frequencies 
specified by the designer are output in crder to aid the 
designer in Baking any modifications to the form of the 
ccupensator that nay be necessary in order to better satisfy 
the system specifications. 

As stated in the previous paragraph six separate cost 
furctions have been incorporated into the function 
sub-program EE (X) . The first three of these are cf the 
error squared type and the latter three are cf the absolute 
error type. These six cost functions represent various 
ccmbinat icns cf one minus the actual response divided by the 
desired response at the discrete frequencies under 
consideraticr , These pre-defined cost functions are 
provided as a matter of convenience in using the program and 
certainly dc net prevent the user from defining his own cost 
function provided this is placed within the FE(X) function 
sub-program. The first cost function combines the deviation 
of both the desired magnitude and phase responses frem the 
actual magnitude and phase responses of the open loop 
system. That is the first cost function is cf the form 



where AKAG and APHASE are the actual magnitude and phase 
response cf the epen loop system and DtiAG and DPHASE are the 
desired iragritude and phase respectively. These quantities 
are evaluated in the cost function and summed ever the 
discrete frequency values in the specified range. The 
second cost function which the user may chocse, uses solely 
the deviation of the actual magnitude response desired 
magnitude response evaluated at the discrete frequency 
values. Thus this cost function is of the general fern 
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AMAG (*) 



DMAG (u.) 



If a minimum phase system is being considered the magnitude 
and phase plots as a function of frequency are uniquely 
related by the Bode theorems [3]. That is, determining the 
magnitude ratio of a system over a given frequency range 
will also result in a unique phase relationship for the 
system. Thus if the system design being considered is of 
the minimum phase type the use of the type two cost function 
described should result in both the magnitude and phase 
specifications being satisfied. The third pre-defined cost 
function uses the deviation of the actual open loop phase 
response fcm the desired phase response in a function of 
the form 



While for ninimum phase systems the magnitude and phase 
responses are uniguely related, it should be noted that the 
fixed less or gain of the system cannot be determined from 
the phase characteristics alone. Conseguently the use of 
the type three cost function may result in a constant gain 
error for the magnitude response of the system. This error 
however is readily observable from the output produced by 
the program and cprrection of this may be accomplished by 
adjustment of the compensator gain value accordingly. The 
ether three pre-defined cost functions are of the same 
general ferm as these just presented with the exception that 
the squared differences are replaced by the absolute value 
of the differences. The general rationale in selecting this 
latter type of cost function was that on occasion the 
sguared type cost function will have a tendency to flatten 
near the Dinimum. With the absolute value cost functicn the 
user may avoid the excessive time required in finding the 
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minimum values when the case arises where the squared cost 
furcticn becomes flat. 

In utili2ing this design program the user should alsc be 
aware that the need for frequency scaling of the problem may 
arise. This is particularly true when high order systems 
involving high frequencies are being considered. The need 
for freguency scaling is primarily a result of the finite 
word size available in the computer system. It must be kept 
in mind that the magnitude response of the system is 
computed as the ratio of two complex numbers and while the 
resulting magnitude may be relatively small the value cf the 
magnitudes of the individual complex numbers forming this 
ratio ma? he large. The situation is further complicated by 
the fact that in computing the magnitudes of the complex 
numbers involved the real and imaginary parts must be 
sguared, thus effectively reducing the maximum power of the 
numbers which can be allowed by a factor of two. If the 
magnitudes of the resulting numbers become too large an 
error message will be printed recommending that the problem 
be freguercy scaled. 



B. DATA INfCT AN l OUTPUT 



In crdei to simplify use cf the program, an attempt has 
been made tc maintain the input information necessary for 
execution cf the algorithm to a minimum, and at the same 
time preserve a certain amount of flexibility within the 
program which may be controlled by the user. In so doing 
seme of the program control variables for various options 
have been set irternal to the program and may not be 
manipulated by the user in the input data deck. These 
internally set variables do not detract from the general 
functioning cf the program, but mention of this is made to 
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bring it tc the user's attention in the event that 
specialized problems might more effectively be investigated 
by resetting seme of these normally preset program 

variables. fiany of the preset items are concerned with 
control of the minimization routine. Such items as the 

upper and lower bounds of the search area, the integer 
programming option, and the seed for the random number 
generator are set to specific values within the main program 
just prior tc calling the minimization subroutine. For 
example for the minimum phase case, the lower and upper 
bounds on the search areas of the compensator numerator and 
denominator polynomial coefficients are set at O.C and 
6 

1.6x10 respectively. For the non-minimum phase case the 

bounds on the search areas of the numerator polynomial 

6 6 

coefficients only are extended to -1.6x10 and 1.6x10 and 

the implicit constraints are adjusted to allcw right half 
plane roots for the numerator polynomial. The integer 
programming option fpr BOXPLX is also preset to a value of 
zero indicating that values for the polynomial coefficients 
need not be integers. In general, selection of the integer 
programming option results in considerably more execution 
time for the program and for the purpose of varying the 
compensator parameters the selection of this option provided 
nothing tc enhance the solution of problems. 

The input data required for program execution and 
control may be divided into general areas as follows: 1.) a 

title card; 2.) a series of cards giving the gain and 
coefficients of the plant transfer function; 3.) a series 
of cards giving the assumed gain and coefficient values of 
the compensator transfer function; 4.) a control card 
describing the plots desired, the type of cost function to 
be used, whether the system contains a zero order hold, and 
the number and rang.e of the discrete frequency points being 
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considered; 5.) cards containing the discrete freguency 
values / the desired magnitude, and the desired phase; 6.) 
finally a card describing the number of iterative trials to 
be allowed in the minimization and whether any output is 
desired from the minimization routine in order to keep track 
of the nanner in which the coefficients are being varied. 
The specific input formats of the above necessary data are 
described ir detail in the following paragraphs. These 
paragraphs are numbered so as to indicate into which of the 
six general areas a particular data card belongs. This 
description coupled with the program listing given in 
Appendix E should provide sufficient information for setting 
up the data deck necessary for execution of the program. 
Figure I V-2 illustrates the data deck used in executing the 
first example problem presented in section C of this 
chapter. In this figure each line of data represents an 
individual data card- The column numbers for the data cards 
are also shewn above and below the data .deck section of the 
figure. 

Title card. 



This 


must be 


the 


first card 


of the 


data deck and 


provides 


a means 


pf 


problem identification if several 


separate 


problems 


are 


run under 


one jot 


name. Columns 1 


through 


46 may be 


used 


to input 


whatever 


alpha-numeric 



characters the user desires for identifying any particular 
problem. This title will also be printed cn the graphical 
output reguested by the user. 

iib Fla ct cain 

The second card contains the gain of the plant in 
columns 1 through 10. If the gain is not explicitly present 
in the plant transfer function then the value input cn this 
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card should te unity. The format used to read the plant 
gain is at E10.0 type format- Thus, the gain may be input 
in either floating point or scientific notation type 
formats. 

2E_. Input form and or der of the plan t transfer function 
numerator,. 



This card specifies whether the numerator of the plant 
transfer function is to be read in factored cr polynomial 
form and what the order of the numerator is. A P or an F in 
column 1 indicates whether the plant numerator data is to be 
read in polynomial or factored form respectively. Columns 2 
and 3 contain the order of the numerator. Thus an input of 
the form EC3 starting in column 1 indicates that the 
polynomial is of third order and is to be read in factored 
form. 

2C_. Plart transfer function n um er ato r. 

The following card(s) contain either the polynomial 
coefficients cr the factors of the plant transfer function 
numerator, depending on which form has been chcsen for 
input. If the polynomial form of input has been selected 
then the polynomial coefficients are input in ascending 
powers of s using an 8E10.0 format. If more than eight 
coefficients are needed they are simply continued on 
successive cards. The coefficient of the highest order term 
is assumed to be normalized to unity and is not read in 
explicity. Thus the number of coefficients actually read 
corresponds tc the order of the numerator polynomial. If 
the factored form of input is chosen then each factor is 
read in cn a separate card with columns 1 through 10 being 
used for the real part of the factor and columns 11 through 
20 for the iaaginary part. Again the format used is E10.Q, 
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thus allowing the factors to be input in either floating 
point or scientific notation form. Since complex factors 
will occur in conjugate pairs the values of these factors 
need only be input once. Thus the number of cards used will 
vary depending upon how many purely real and how many 
complex factors are involved. 

2E_. Input form and order of th e plant transfer function 
denominator . 



Again the form of input is specified as in card number 
2E along with the order of the plant transfer function 
denominator. 

2E_. Plant tra nsfe r f unct ion de no minato r. 

The format used to input the denominator is identical to 
that used fcr input of the numerator contingent upon the 
desired form specified for the input. That is to say, it is 
possible to input the numerator in polynomial form and the 
denominator in factored form, or vice versa, since each time 
a polynomial is input as data the input form is uniguely 
specified . 

li-s Assumed compensator tr ansfer function g ain. 

Following the input of the information about the plant 
the description of the form of the compensator must be 
supplied. This is begun by inputing an assumed compensator 
gain using the same format as was used to specify the plant 
gain. Cue to the manner in which the compensator transfer 
function nuierator coefficients are handled within the 
program, if the assumed compensator gain is set to zero, for 
lack of any better value, whis will result in the starting 
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guess of all the numerator coefficients being zero. This 
will not cause the program to terminate prematurely, but the 
initial guess of zero for all the coefficients may not be 
what is desired by the user. If doubt exists as to what 
value gain should he assumed, a value of 1 is suggested as a 
convinient arbitrary starting point. 

3 B_. in put fcrm and order of the assumed ccm pensa ter 
nun era to r . 

The netted for specifying the form of the input fer the 
compensator initial numerator parameters is identical as 
that used in specifying the plant numerator. A P indicates 
polynomial fcrm and an F indicates factored form. 

3C_. Assumed initial values of c ompensa tor n umera tor 
para mete rs. 

In order to maintain continuity in the input format the 
scheme and format used to input the compensator numerator 
parameters is identical to that used in specifying the plant 
transfer functior numerator. 

3C_. Input form and order of the assumed cc mpensa ter 
de rernina ter . 

As in previous cases either polynomial cr factored form 
is specified along with the order of the compersator 
de remina ter . 

31 . Assumed init ial values of comp ensato r d enomi nator 
parameters. 



Once again the same format as used to input the plant 
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data is used in prder to maintain the continuity cf the 
input foraatirg. 

4 Pro era ii con trol car d. 

The program control card is used to input various 
quantities controlling both the range and number of discrete 
points used in the minimization of the cost function, the 
output desired from the program, and also whether a zero 
order held is present in the error channel of the system. 
As many as twelve separate quantities may be input on this 
single control card in order to direct execution cf the 
program. Columns 1 through 10 specify the minimum freguency 
(WMIN) , in radians, of the total range being considered. 
Columns 11 through 20 specify the maximum value (WMAX) , in 
radians, cf the range of frequencies being considered. Eoth 
of these values are read using a 2E10.0 format. Following 
the maximum and minimum freguency values, up to ten input 
quantities are read using a 1013 format. Columns 21 through 
23 specify the number of discrete frequency points ( NOMEG) 
over the given range at which the open loop freguency 
response is to be determined and the cost furction 
evaluated. The next input quantity (KNOW) , in columns 24 
through 26 specifies whether the discrete frequencies of 
interest are to be read from data cards or automatically 
computed. Normally this quantity is set to unity indicating 
that the discrete frequency points of interest will be read 
from data cards. A value of zero indicates that the 
discrete freguency values are to be incremented linearly by 
the program starting at the minimum freguency value. A 
value of two indicates that the discrete frequency points 
will be incremented logarithmically begining with the 
minimum freguency value. Columns 27 through 29 indicate 
whether a Ecde plot of the results is desired. If this 
variable (NECEE) is set equal to one no Bode plot will be 
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output. A value of zero for N BODE will result in fcoth 
magnitude arc phase plots. Columns 30 through 32 control 
the output cf the Nyguist plot. If NYQS1 is set to unity 
thee the Nyguist plot will he suppressed. A value cf zero 
for NYQSX will result in a Nyguist plot output, however this 
is plotted crly over the range of freguencies between WMIN 
and UMAX. Columns 33 through 35 indicate the presence of a 
zero order hold (IZOU) in the error channel of the unity 
feedback system. .With IZOH set to zero the system is 
treated as being of a continuous nature. With IZOE set 
egual to one a zero order hold is considered to be present 
in the error channel and the open loop magnitude and phase 
computations are modified accordingly. Columns 36 through 
38 (NICHCI) specify whether a magnitude versus phase plot of 
the results are desired. As with the other plot options, a 
value of zero will result in a plot being produced while a 
value of ere will suppress the magnitude versus phase plot 
as part cf the output. Columns 39 through 41 contain the 
variable indicating whether non-minimum phase solutions are 
to be considered. With this variable (NHINFS) set egual to 
zero only minimum phase solutions will be allowed. A value 
of one indicates that a non-minimum phase solution will be 
allowed provided that it is the one which minimizes the cost 
function. Columns 4 2 through 44 indicate if the desired 
magnitude profile values to be input are in decibels. If 
this variable (IDB) is zero the magnitude values to be read 
are dimensionless gain values, while a value cf one 
indicates that the magnitude values to be input will he in 
decibels. Columns 45 through 47 contain the variable IPLCT. 
This guantity is used to choose between printer plots or 
Calccmp plots for the graphical output. A value cf zero 
will result in printer plots and a value of unity will 
provide Calccmp output of the graphs desired. It is 
recommended that the printer plot option be used in 
initially selling a particular problem since this type of 
output requires less turn around time than the Calcomp 
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output. Finally columns 48 through 50 specify the type of 
cost function (ICOST) the user wishes to choose among the 
six preprogrammed functions available in subroutine FE (X) . 
The value specified should be between 1 and 6 corresponding 
to the particular cost function desired. If this value is 
zero the program will default to the type one cost function. 

4E_. Sa mpling period^ 

This card must be inserted only if a zero order held is 
specified in the error signal channel. The sampling period 
T, in seconds, is specified in columns 1 through 10. If the 
system being considered is of the continuous type then this 
card must be emitted. 

5JU Di screte frequency values consider ed . 

Next fcllcws a set of data cards specifying the values 
in radians per second of the frequencies at which the 
desired gain and phase values of the open loop freguency 
response are specified and at which the cost function is to 
be calculated. These values are input in an 8E1Q.0 format 
using as many cards as necessary to specify all the 
freguency values at which the program is to execute the 
algorithm . 



5E_. Gain profile of the desired open loop response. 

This set cf data cards contains the desired open loop 
gain response at the discrete frequency values previously 
specified. The same 8E10.0 format is used here as was the 
case in reading the discrete frequency values. The same 
number cf cards should be used as was used in specifyirg the 
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discrete frequency values. 

JC.. Phase pr ofile of the desired r espon se . 

Here fellows a series of cards specifying the desired 
open locp phase response, in degrees, at the discrete 
freguency values being considered. 

JL Minis izaticn trials card. 



Here both the number of trials to be allowed in the 
minimization process and the print interval for diagnostic 
purposes are specified. These quantities are input using a 
215 format. If the number of trials is zero then a default 
value of 2CCC is assumed. If the print interval is 
specified as zero then no output from the minimization 
routine will result. Caution should be exercised in 
specifying the print interval in that if this interval is 
made too snail an excessive amount of printed output may 
result. The program user is referred to the comment section 
of subroutine BOXPLX in Appendix A for guidelines in 
selection of the print interval. 
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Figure IV-2 



c. 



SERIES COMPENSATOR DESIGN EXAMPLE PROEIEMS 



In this section a number of example problems are 
presented and discussed in order to illustrate the use cf 
the computer aided compensator design algorithm- In the 
hcpe cf presenting a clear, concise illustration of the 
algorithm and its use the example problems begin with 
relatively simple straight forward text bock examples and 
gradually progress to problems of a slightly more difficult 
nature- Also included are problems that are intended to 
illustrate some limitations of the program in that these 
problems were not completely solved by the algorithm cr they 
may have reguired a certain amcurt of manipulation cf the 
input data in order to obtain a solution- These latter 
types of problems have been included since it is important 
that the designer using any computer aided design algorithm 
understand the algorithm’s limitations and peculiarities as 
well as its intended capabilities . It is net intended to 
mislead the reader by implying that the example problems 
presented here provide a comprehensive presentation of all 
possible limitations of this automated design program. 
Rather they represent those limitations which were 
encountered ir the process of verifying the program ever a 
necessarily limited period of time- With this in mind the 
first example problem which deals with a relatively simple 
lead compensator is presented. 



1- Compe nsa tor De sig n , Ex ample _K 



Consider the uncompensated system shown in figure 

IV-3. 
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Figure IV-3 



The asymptotic Bede diagram for the open loop system is 
shown in ir figure I.V-3A. As can he seen from the open loop 
Bode plot / ever theu-gh the system is stable it is lightly 
camped. A single section lead compensator is to be used to 
increase the phase margin. A complete discussion of the 
reasons for selection of the lead compensator and a detailed 
discussion of this particular example is presented by Thaler 
and Erown [18]. The transfer function form of the single 
section lead compensator consists of a first order numerator 
ana a first order denominator. A desired magnitude and 
phase profile, consisting of ten discrete frequency points 
over the range from 0.2 to 20.0 radians, was selected to 
give a phase iiargin pf approximately 45 degrees in order to 
increase the damping of the system. As an initial guess the 
ccnpensator was assumed to have a gain of 1.0, a zero at the 
origin arc a pole at -1.0. This information was supplied to 
the program in the prescribed format and the optimized 
results in the form of magnitude and phase plots are shown 
in figure IV-3E. Here, as throughout the remainder of the 
example problems, the diamond shaped symbols represent the 
desired values specified by the program user, the X symbols 
represent the values computed by the program at the discrete 
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frequency values specified using the compensator parameters 
returned as the solution from the minimization routine, and 
the continuous curve represents the response computed by the 
program again using the parameters returned as a solution 
but with a mere dense set of freguency values used in order 
to show any resonant peaks or anomalies in the magnitude and 
phase curves which might be difficult to visualize from a 
sparsely selected set of discrete values. As can be seen in 
figure IU-3E the results match the desired response well 
over the range of interest. Figure IV-3C shows the same 
results using the magnitude versus phase plot to portray the 
open loop frequency response. The numerical values of the 
compensator parameters returned from the program as a 
solution along with the starting values used are shewn in 
figures IU-3E, IV-3E, and IV-3F. A comparison with the 
solution values given in ref. 18 will show that they agree. 
Figure IV-3 G shows plots of the differences between the 
actual ard desired open loop frequency response of the 
resulting system at the discrete frequencies considered in 
the minimization process. 

This problem as described required approximately 500 
iterations within the minimization routine to achieve the 
solution shown. This however includes two restarts to 
attempt insuring that a global minimum had been achieved and 
that termination was not due tc a local minimum. From the 
begining tc the first restart required approximately 250 
trials with the cost function being reduced from a value of 
6.23 to a value of 0.00711 (recall that the minimum 
achievable under ideal condition is zero which due to 
machine errors alone in representing numbers would not 
likely ever be achieved). The entire solution required 32.5 
seconds of CEO time, a good 10 seconds of which was used in 
the production of Calcomp plots. These numbers concerning 
the time required for solution are presented for what they 
are worth in that extrapolation of these numbers tc more 



51 



explicated problems is not a straight forward task. As 
mentioned in Cbapt III the minimization routine does become 
inefficient when the number of variable parameters is 
allowed tc become large. It was felt that a detailed 
investigation of the efficiency and time required for 
solution cf problems as a function of the number of variable 
parameters and the number of discrete frequency values tc be 
considered was beyond the scope of the basic developmert of 
the program. However, the amount of time required for 
solution cf seme of the example problems that follow will be 
given in order to give the user an "intuitive feel" for the 
approximate amount of time that may be need to handle 
certain types cf problems. 

In order to illustrate that the starting point 
chosen should have reasonable values, but that the exact 
values chosen should not effect the final results 
signif icartly , the same problem was run a second time with 
the pole and zero of the compensator both chosen tc be at 
the origin. Again an assumed initial gain cf one was chosen 
for the compensator. The graphical results cf the solution 
are shown in figures I7-3H through IV-3J. The numerical 
results returned using these new starting values are not 
shewn here since the.y were identical to within three decimal 
places with those values shown previously. 
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Uncompensated Open Loop System Bode Plot 
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Magnitude vs . Phase 



Figure IV- 3 C 
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Computer Numerical Output Example 
Figure IV-3D 



■[ o 1 00* T - = NMV'J NCiiJ. IlNfH >f-.i J V-fJV >U kMO-> *V-J 7 Id I IcM 



! 



w_# 

-— rvi r^i rv i -\j rs, rv r\i rs; .Ni 

<rcccccoccao 
X . 



O.G; i U’U--l-'Li *: Li 'ij' ij 




c'oooococoo 


f | 


□ c'coroocooo 




J'CCOSCCOO o 




Ol j* C O C DOCOO C 




— CCi— J' 




v) o OCrviro^r-orsir*' 






1 1 


C O' i^<— >r\jr j *n) 

l l 1 l l 1 i I i i 

i 


1 1 
1 1 


i ; 

rvi Gvi rvi r\j f aj cn; 


1 1 


i ccorooooo: 


1 t 


7T 


1 t 


^ ro li"' 


1 — 1 


G^OJ^T'GvNNG 


1 LU | 


. ■ ” gg cr- c 


1 H- 1 


i L •*- <- -C ;5 — O ' 




(.0 o — CG c 


* 1 


< OOO <V !"**, U“ O 'H*”"* 


j r ; rH 






C- r — — «— 


! ^ 1 ® 


i i i i i i i i i i 






- & 


"■*" 


1 g 


O 






• - i 3 


-OCCCGCGOCG 
O I 


! ^ ! m 


.. '.c ^ 


l T • 4 J 


c; ■ H Cr. f\: «\J r~ o >J* ^ j j 


1 O 1 c pp 


“ v-r^LT. .'Hs-Cir^ 


• — i jsum 


, — . _» i~~ rv C «T* 'V 


1 1 4 -> 1 


c o <* • '* 


! - 1 3 > 


O' C O O C* 


1 ^ 'Oh 


-r.'.C cr- 






• ^ ' rH <D 


, . — 1 *— — ' p. -'G ? — “ C "•"* *4* 


1 r i <3 u 


d 1 l 1 l 1 


r": o p 




1 1 -r| {J! 




^ 1 5-1 -H 




' G ! 3) ^ 





« <■ 1 g 


c O OO O o GCOO' 




O 1 


\ < \ ^ 


• . 


1 X 1 


ii. r.— — • — * _? v ■ r *• *\i c z* 


■ n- » M 


ZL “ O i' ~G vT -sT ’■*” ~T> T 


i I (D 


C “'.1 C vJ ~ — * Vi O — • 


» ~ » +J 


h g • ■ . f— •'••; 


i X. . 3 


•“v/'vTi^' v' C '7 .j n — * CG 


i ) 1 • Q, 




* ~ i g 




1 - -L , 0 


^ r— « ~ J" 


» . u 


. i • * t ■ 


i ^ » 


i ■ i 


1 v- 1 


j ; | 1 


1 Vi'-O 1 




1 } — t 


O C/ o o o o o O o c 

till 


1 X 1 

ll-il 


* v*M 


1 1 — * — 1 


>Q?OO vC o o o o 


1 OV'I 1 


! GOCCO-O'J'OOC 


1 'J>- • 


— O •_ O 3 j“- O ^C- 


1 1 


— 6 - r ^3 cgu 5 o c 


1 Ui..J 1 


C O r"*- _-i ' «j IT x-O — « .Vi c 


1 X X i 




1 1 — I 


Cl. v: 'C O —• 'NJ -V ^ — • VI 


1 1 




1 r~ l 




! (— ! 



58 




o 

o 




- 








- 




A 
l — i 


A A 


.02° 




A 










A i 

1 1 1 1 1 III 1 


1 1 1 1 1 


, frequency 
Mil i in 


°°0. 


1 


1.0 




(rads/sec.) 










A 


.02° 




A 






O 

o 




Phase Difference vs. 

Figure IV-3G 


Frequency 

A 





59 





60 




Example 1, Run#2 
Figure IV-3I 
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Compensator Ees iq n, Example 2 a 

Ibis example problem (discussed in detail in ref. 10) 
has been chosen net only to familiarize the reader with the 
use of the program, but also to illustrate one of the 
shortcomings of the design program. In particular, while 
the program was successful in finding the correct solution, 
difficulties were encountered in the root finding 
subroutines. In this particular problem, during the final 
phases of execution the roots cf a fourth erder polynomial 
must be feure in order to calculate the phase response using 
subroutine SlbULA which does these calculations using a more 
densely spaced set cf discrete frequencies than it the 
minimizaticn routine. Since numerical root finding routines 
have difficulty in finding roots for fourth order 
polynomials, the lack of good convergence for the roots of 
the system in both rcot finding subroutines caused the phase 
calculations normally done by subroutine SIMULA to be 
aborted. It should be emphasized that failure of the root 
finding subroutines in this case only effects the generation 
of the "continuous" phase response and in no way created any 
problems for the minimization part of the program as 
evidenced by the resulting correct values for the 
compensator parameters. The original uncompensated system 
being considered is shown in figure IV-4. A Bode diagram of 
the uncompensated system (figure IV-4A) reveals that the 
system as it exists is unstable. If stability were the only 
consideration a simple attenuation of the plant gain would 
be sufficient to provide adequate phase margin to insure 
this. However, as with many series compensation problems, 
here we must effectively strike a compromise between 
effectiveness cf control (generally associated with a higher 
gain) and stability pf the system (generally associated with 
a lower gain). As discussed in ref. 10 a simple gain 
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attenuation wculd result in stability being achieved, but 
also has at undesirable effect upon the velocity error 
constant. Thus as discussed in detail in ref. 10, a lag 
compensation retwcrk appears as the best candidate to 
accomplish tfce job. 




Figure IV-4 

In this particular problem twenty discrete frecuency 
values over the range from 0.005 to 50.0 radians were 
selected ara the desired magnitude and phase values 
specified at these frequencies. As in the previous example 
problem the form cf the compensator is selected as having a 
first order numerator and first order denominator, that is 
in this case a sincle section lag compensator is to he used. 
A type one form of the cost function was selected for the 
minimization process. The magnitude and phase profiles of 
the desired response along with the resulting response with 
the compensator included in the open loop system are shown 
in figures I\-4E, IV-4C, and IV-4D. As mentioned, the 
continuous curve for the phase response was not computed 
because cf the failure of the root finding routines to find 
the open lccp system roots. It can be seen in figure IV-4C 
that the phase response does agree well at the discrete 
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frequencies specified. Since we are dealing with a minimum 
phase system and figure IV-4B shows that there exists no 
resonant peaks in the magnitude curve the user may conclude 
that the jhase response will also te well tehaved fcr the 
various freguency values between the explicitly specified 
discrete values plotted. This case should serve to 
illustrate that the program reguires that the user still 
provide ar interpretation of the results. The computer 
output cf figures IV-4E, IV-4F, and IV-4G shows that the 
compensator parametexs are in agreement with the problem 
solution as presented in ref. 10 if the reader wishes to 
make the comparison. Figures Iv-4H and IV-4I show plots of 
the differences between the actual and desired magnitude and 
phase at the discrete frequencies selected for the 
minimization process. 
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Figure IV-4C 




Figure IV- 4 D 
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Computer Numerical Output Example 
Figure IV- 4 F 
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Phase Difference vs. Frequency 



3. 



Compensator De sig n , Example 3_. 

Another lag type compensation problem is presented 
in this example. This time, however, the existence of a 
fourth order polynomial for either the open loop system 
transfer function numerator or denominator has purposely 
teen avoided tc insure that the program will run tc 
completion. The uncompensated system block diagram is shown 
in figure 17-5. 




figure IV-5 



A Bode diagram for the uncompensated system along with a 
detailed discussion of the problem may be found in ref. 17. 
The uncompensated system does not possess an adequate phase 
margin tc insure proper system performance. Thus, a single 
section lag compensator is to be used in order to achieve 
the desired open loop frequency response. Desired gain and 
phase prefiles are selected for the compensated open loop 
frequency response using 20 discrete freguency points over 
the range frem 0.01 to 100.0 radians. The prefiles 
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selected along with the results achieved are shewn in 
figures I?-52, IV-5B, and IV-5C. The numerical values 
assumed as initial estimates of the compensator parameters 
and the numerical results of the output are shown in figures 
IV-5D, IV-5I, and IV-5F. Piets of the differences between 
the actual magnitude and phase, and the desired magnitude 
and phase are shown in figures IV-5G and IV-5H. In this 
particular prcfclem approximately 1200 iterations of the 
minimization routine and 1 minute and 25 seconds were 
reguired for solution of the problem. In order to 
illustrate the idea that for minimum phase systems the 
compensator parameters may be determined to within orly a 
fixed loss or gain when using only the phase profile, this 
same problem was solved again using the type 3 cost furction 
which considers only the difference between the desired and 
actual phase at the discrete frequencies. The results are 
shewn in figures IV-5I, IV-5J, and IV-5K. As can be seen 
from the magnitude plot in figure IV-5I, there is a constant 
gain errer in the resulting magnitude curve. The numerical 
values assumed at the start and those returned after the 
minimization cf the cost function are shown in the computer 
output of figures IV-5L, IV-5M, and IV-5N. The resultant 
magnitude response is higher than desired. The necessary 
correction that must be applied to obtain the desired 
response may be read directly from figure IV-50, which shows 
the difference between the actual and desired magnitude 
response. That is, as can be seen from figure IV-5C, the 
resultant magnitude of the open loop system after 
compensation is approximately 36.8 db higher than what is 
desired. Since the program does not alter any of the plant 
parameters this means that the compensator gain value 
returned from the prpgram when the type 3 cost function was 
used must be decreased by approximately a factor of 69.2. A 
comparison of the compensator gain values shewn in figures 
IV-5I and IV-5N will show that they due indeed differ by 
this value. Also the astute observer may notice that there 
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is a slight difference in the compensator pole value 
returned as the solution when the problem was solved using 
the type 3 cost function. While the difference in values is 
small it is suspected that in this particular case it is due 
to the fact that only a finite frequency range has been 
considered acd the lower frequencies have been "slighted” 
somewhat in that the nonlinear phase curve has net yet 
flattened cut at 0.1 radians. This is also the same 
situation that existed the first time the program was 
executed, but in this case the magnitude prefile was also 
taken into account in the cost function and any deviations 
of the resultant magnitude values dominated the cost 
function value. In other words, it appears that in this 
particular problem the cost function is more sensitive to 
variations in the magnitude values than in the variatiens in 
the phase values. 
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Figure iV-5B 




Figure IV-5C 
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Computer Numerical Output Example 3 , Run #1 
Figure IV-5D 
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Magnitude Difference vs. Frequency 
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Phase Difference vs. Frequency 
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Figure IV-5I 




87 




88 



Cvl 

o 



1 


o 


1 








O 




1 


UJ 


• 


1 










1 


c 




1 


j 




C 






0 


1 








1 




1 


c 






1 




O' 




1 


0 


i 


i 


VO 




i 




oo l 




j 




1 








1 


0 






l- ! 




1 




h- i 


6 


1 




CtZ 




QC 


* 


1 


0 


| 


i 


O i 




1 




CL I 


• 






0 


0 


a 


i 


1 


0 




1 


O 




• 




O 1 


W-t 


Cl 




K 


O 


0- 




1 


0 


; 




CL 




CL 




CM 1 




o 




< 




<3 


1 


1 


• 










<-7 




1 




►— • 




21 


UJ 


2 






rH 


CL 




CL 




£■ 




1 


II 


< 




y—* 


0 


>— ■ 




1 




O 




O 




< 




PQ l 




cl 




21 


0 


2; 




1 




h- 




1“ 




2 




CO 1 


2: 


UJ 




CO 


O 


Oh 




1 


11 


<r 




< y~ 




>_« 




1 


>— i 


y ' 




2 


© 


20C 




1 




a 




Cl CL 




y 




LU 1 


< 


h 




UJ 


s-> 


UJ< 




1 


2 


UJ 




UJ< 




O 




_J 1 


o 


2 




Q 


6 


QO. 




1 


y—t 


2: 




SI Cl 




2 




Cl 1 




oo 




VO 


• 






1 


<t 


©VO 




2 




lUOO 


j 


s: i 


2: 


2 




2 


■—4 


2 > 




1 


O 


2 


1 


Z> 




Q 




< i 


o 


OLL 




GLL 




OQC 




1 




LL 




CL 




LL 




X | 




~o 




— O 




*-< 




1 


2 


20 




2<r 




20 




LU I 


\- 


K 




►— 




1— 2 




1 


O 


C 




02 




O 


j 


1 


o 


LJOO 




OOO 




O — 




l 


►— * 


‘—'O') 




►— < •— « 




— O'! 




• | 


2 


20C 




2 <X 




2 0 


0 


© 1 


h- 


►— LL 


0 


1— O 


O 


K- O' 


O 


M 1 


13 


Z>oj 




Z)LU 


O 




• 




O 


OUJ 


0 


O <t 


• 


O'JJ 


O 




LL 


U.-S 




LL .2 




LL 5 


0 


© 1 


2 


2 « 




2 >_ 


rj 


2 3 




5 L 1 




LJ 




O 


UJ 


►— « 




1 


O 


-JO 


UJ 


_>•“« 




_>o 


LU 


M | 


<x- 


CCLL 




QC LL 


O 


a 






U. 


LL Ok. 


0 


lL 




LL CL 


O 


r— I 


UJ. 


LU 


J 


UJ 


O 


UJ 


1 


: • 






w 








LJ 


CL 1 


LL: 


LUCJ 




LL O 


O 


LL 


! 


i 1 


QL 


iXO 


0 


a ] 




ao 


O 


O 1 


00 


00 2: 


i 


O') 2 


; a 


VO 


i 


j I 


LU 


LU 2 


i 0 


UJ ! 




LU2 


CJ 


1 


l; 


2.M 


i 




LU 


2 


i 




U. 


u_ *— « 


! 0 


LL. 1 


i 


u. >— 


w 


• 1 




<0.0 




<a LJ 


O 


<t 


1 

: 




VI 


1/1 OJ 


0 


oO 




VJLJ 


0 


2' 1 


cc 


a: 2 




Cl 2 


• 




] 


© 1 


2 


22 


j J 


2 


O 


22 


• 


LU | 


i— 


h-uj 


I* 


^-*LU 




>- 


1 


<j l 


< 


•-f UJ 






O 


<1 LU 


•—1 


CL 1 


; 


O 


i 


O 




1 L— 


! 




CL 


QCO 


J 


X I- 1 


I 


QCO 




X 1 


C 


Q0O 


1 


CD OO 




Ooc 


j 


UJ | 


H- 


K.LO 


j 


1 — CL \u 


►-TOO 




U 1 


LU 


iu <1 


1 


UJ <1 




IU<J 


1 


0 1 


1 


<1 


j 


< 


O 


' <1 




LJ 1 


H- 


>— ! 


! 


H- 




f— <-L 




0 1 


(X 


IX 


! 


oc a. 


O 


oc 




1 


<l' 


<12 


O 






■d 


i 


© 1 


0 


G2 


•0 


O 


O 


02 


0 


1 1 


O 0 


OOt— « 


CJ 


oo»-« 




V)J 




CJ 1 


►— 


►— *— * 


0 




O 


(->-• 


0 


1 1 


2 


2 




2' 








w 1 


< 


<l 




<1 <i 


O 


<l 




1 I 


UJ 


LL .00 


UJ 


UJ 




LU UJ 


rj 


LJ 1 


00 


coijn 


1 III 


*/UU 1 


O 


00 00 


UJ 


1 


a 


CLh- 


CJ 


CL ►— 




a. oc 


• 


• | 


2 


2H- 


© 


2 CL 




2 ►— 


0 


UJ | 


20 


XZ 


LU 


2l2 




2 . 


CJ 




UJ 


UJ2 


LJ 


LU 


—4 


UJ 2 


, 0 


_J | 


O 


Ouj 


O 


Ouj 




O 




1 1 


CL 


CLUJ 


O 


CL 


| 


CLUJ 


0 


h- 1 


o 


o>-* 


O 


LJ*-* 




0 






2 


2.*-* 


O 


2. 




2l*-* 


0 


M 1 


2 : 


zu 


O 


20 




2 






CO 


OO 


O 


CJ 




OO 


0 


1— 1 


ZD 


JM 


LJ 


_>•— 


O 


-J 




1 


0 


O — 


O 


O 




O— ' 


CJ 


1 




u_ 


• 


LL 


• 


• • 


1 






LL 


• 


«• 




lL 


• 


1 




u_ 


f-4 


LL 


O 


UJ 




| 




*-L 


“1 


LU 




LL 


<— 4 


1 




LU 




UJ 




ac 








UJ 




iX , 




lL 




1 




LJ 




rj 




<1 




1 




0 




<1 




O 




1 

1 


1 


O 




0 








1 




0 








O 





! 

O 

o 



89 



Computer Numerical Output Example 3, Run #2 
Figure IV-5L 
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Magnitude Difference vs. Frequency 
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Phase Difference vs. Frequency 



4. 



Cc a fensa t cr Desig n , Example 4 

In this exanple a slightly more complicated problem is 
presented, in that the final solution reguires a double 
section compensator of the form of a notch filter. A block 
diagram of the uncompensated system is shown in figure IV-6. 
As discussed by Thaler and Erown [18], from whom this 
example was taken, the uncompensated system is unstable and 
the Ecde diagram of the uncompensated system shown in figure 
IV-6A indicates a phase margin of approximately -30 decrees. 
With such a large negative phase margin for the 
uncompensated system and the rapid drop in slope of both the 
gain and phase curves the designer might suspect that 
compensation using only a single section compensator will be 
difficult to accomplish. Again desired magnitude and phase 
profiles are selected and initially a single section of 
compensation is assumed to observe how close to the desired 
response the system will perform. As can be seen in the 
resultant graphical putput of figures IV-6A, IV-6B, and 
IV-6C, the single section compensator fails to meet the 
reguired frequency specifications. In fact, while it is 
possible to stabilize ths system with a single section 
compensator the bandwidth will be excessively large and the 
desired magnitude profile prevents the program from 
accomplishing this. The numerical results returned frcm the 
single section compensator run are shown in figures IV-6D, 
IV-6E, and IV-6F. As can be seen in figure IV-6F, the South 
test of the characteristic equation also indicates system 
instability which is to be expected. The magnitude and 
phase difference curves of figures I7-6G and IV-6H, which 
indicate the difference between the specified magnitude and 
phase profiles and the values actually achieved, indicate 
that additional compensation may be needed at the higher 
freguency ranges. Thus a double section compensator was 
assumed and the results achieved are illustrated in figures 



94 



IV-6I, I?-60 # and IV-6K. The program is able to satisfy the 
specif icaticns with a double section of compensation and the 
parameters reguired to accomplish this are shown in figures 
IV-61 , IV-6t, and IV-6N. The differences between the 
desired and specified magnitude and phase values fcr the 
double section cf compensation inserted in the systen are 
shewn in figures IV-60 and IV-6P. 
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Figure IV-6 
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IV- 6 C 
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Computer Numerical Output 
Figure IV-6D 
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Computer Numerical Output 
Figure IV-6E 
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Computer Numerical Output 
Figure IV-6F 
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(rads/sec) 
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Phase Difference vs. Frequency, Single Section Compensator 
Figure IV-6H 
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Magnitude vs. Frequency, Double Section Compensator 

Figure IV- 6 I 




106 



Phase vs. Frequency, Double Section Compensator 

Figure IV-6J 




Magnitude vs . Phase , Double Section Compensator 
Figure IV- 6 K 
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Computer Numerical Output 
Figure IV-6L 
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Computer Numerical Output 
Figure IV-6M 
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Magnitude Difference vs. Frequency, Double Section Compensator 

Figure IV-6Q 
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ase Difference vs . Frequency , Double Section Compensator 
Figure IV-6P 



5. 



Compe nsa tor D esign, Ex am ple 5^ 

lhis final compensator design example represents a 
prcfclem cf a slightly different nature than those presented 
to this point. The servo system under consideration here 
has a zero order held present in the error channel ard in 
addition, due to the large coefficient values that result in 
the plant transfer function, the problem reguires frequency 
scaling in cider to insure that the numerical limitations of 
the particular computer that was being used would net be 
exceeded. The original single loop system is shown in the 
blcck diagram cf figure IV-7. Series compensation of this 
system was desired in order to meet the following 
specifications; 1.) the open loop magnitude at 377 
radians/seccnd > 3£db, 2.) the gain crossover frequency > 
1865 radians/second, and 3.) the phase margin at gain 
crossover > h5 degrees. These specifications were to he met 
using only a series compensator and the parameters of that 
part of the system shown in figure IV-7, other than the 
compensator, was tc remain unchanged. An initial attempt at 
using the CAIICO program to design a compensator resulted in 
a warning of excessive number size during computation and 
recommended frequency scaling of the problem. Thus a scale 
factor of 10C0 was selected as a convenient value and the 
entire system, including the sampling frequency of the zero 
order hold, was scaled down in frequency by the value of 
this scale factor. Using Laplace transform notation, this 
was accomplished by a straight-forward substitution of 
s = 1 000 £ , where S represents the new scaled freguency 
values. £ blcck diagram shewing the values of the system 
parameters after scaling is given in figure IV-7A. When 
this type of scaling is performed the resultant numerical 
values cf the program need only be multiplied ty the 
appropriate constant value in order to be corrected for use 
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in the unsealed system. 



Originally, work on this particular system indicated 
that a fifth order compensator could provide the necessary 
phase margin and gain crossover frequency, tut that the 30 
db specification at 377 radians/second was not satisfied. 
Thus, with this information as a starting point the 
magnitude profile of the system with the fifth order 
compensator in the loop was adjusted at the .low freguency 
end to satisfy the 30 db requirement at the scaled freguency 
value of C.377 radians/second. The phase prefile was left 
unchanged. The results of the program can be seen in 
figures IV-7E through IV-7I. As can be seen from the 
graphical output the desired specifications have been 
satisfied by the resulting compensated system. 
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Unsealed System Block Diagram . 
Figure IV-7 
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Scaled System Block Diagram 
Figure IV-7A 
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Magnitude vs. Frequency 
Figure IV-7B 
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Figure IV-7C 




Magnitude vs . Phase 
Figure IV-7D 
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Computer Numerical Output 
Figure IV-7E 
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Computer Numerical Output 
Figure IV-7F 
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Computer Numerical Output 
Figure IV-7G 
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Magnitude Difference vs. Frequency 
Figure IV-7H 



o 

< 1 ) 

W 

\ 



% 

rd 

U 




124 



Phase Difference vs. Frequency 
Figure IV-7I 



V. 



SYNTHESIS OF TRANSFER FD NCTIO NS FROM FREQUENCY 
RESPONSE DATA 



A. GENERAL DISCUSSION 



Id tfce process of running various test problems to 
exercise and verify the algorithm presented in the preceding 
chapters, it was found that in addition to its use as a 
ccapensatcr design program, the basic method could also be 
applied tc the synthesis of transfer functions from 
freguency response data. That is, if the magnitude and 
phase measurements of some system are obtained it is 
possible tc use the program (CALICO) in order to obtain an 
approximate linear analytical expression for the transfer 
furcticn cf the system from which the measurements were 
obtained. This type of problem is quite often encountered 
in^the modeling process of dynamic systems. That is, the 
engineer may have measured input-output data of the system 
in the form cf a freguency response and is confronted with 
obtaining a transfer function which will model this system. 
This may be the result of the process being toe complex to 
be analyzed and the pertinent equations written on the basis 
of physical laws or the governing physical laws may net be 
completely understood in the case of systems involving new 
technologies or reasearch efforts. Quite often the 
procedure followed in synthesizing or determining an 
approximate linear transfer function from frequency response 
data is tc construct a Bode plot for the measured data and 
to fit the measured characteristic with straight line 
asymptotic approximations where possible and guadratic 
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factors where sharp peaks in the measured data occur. As in 
compensator design by classical frequency domain techniques, 
this procedure for synthesizing transfer functions can 
involve a considerable amount of time consuming and tedious 
"cut and try 11 cn the part of the engineer. 

The algorithm previously presented is so structured that 
the use cf the program for this type of transfer function 
synthesis is possible with no modifications to the coding 
itself. Since the algorithm was not originally developed 
for the specific purpose of transfer function synthesis 
there are several precautions that must be taken in order to 
insure that the transfer function obtained does indeed 
adequately represent the system from which the measurements 
were obtained. In particular when using the program in 
order tc synthesize a transfer function, a separate 
simulation program for time domain response should be used 
to verify the transfer function obtained. This procedure of 
checking the time domain response should be followed fcr any 
transfer function obtained from frequency response data. 
The primary concern here is that it is possible to 
accurately represent the frequency response cf some systems, 
over a finite frequency range, by a transfer function of 
higher order than the true transfer function, but in some 
cases the dominant response to standard time domain inputs 
will be considerably different than that of the actual 
system. Also, while there is a direct relationship between 
the time domain and frequency domain responses through the 
inverse fourier transform integral, this integration, 
strictly speaking, is taken over the entire frequency range 
from -oo tc *oo. i n practical applications it is obviously 
necessary to limit the range of frequencies to be measured 
to some finite region. If, however, this region is not 
sufficiently large for a particular system the resulting 
transfer function obtained from the frequency response data 
may not accurately describe the system. If this is the 
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case, it should become readily apparent in a time domain 
simulation cf the equation obtained from the freguency 
response data. The possible pitfalls in the use of this 
algorithm tc accomplish transfer function synthesis, such as 
the one just mentioned, will be discussed later when some 
example problems will also be presented. 

A linear time invariant system may be expressed as a 
ratio cf two frequency dependent polynomials of the form 

n n- 1 

a s + a s + ... ♦ a s + a N (s) 

n n-1 10 



o m-1 

bs ♦ b s ♦... +bs+b D(s) 

m m- 1 10 

where s is the frequency dependent complex variable cf the 
Laplace transform. Thus the value of this transfer function 
is dependent upon the value of s and the unknown 

coefficients cf the numerator and denominator. Note that 
this is analogous to the situation that occurs in the design 
of a cascade compensator for a unity feedback system with 

the plant represented by a unity gain block. Thus fcr the 

purposes cf using the program, if we imagine the measured 
frequency response data to represent some desired system's 
open loop frequency response profile, and represent the 
plant of that system by a constant gain of unity, then the 
transfer function H(s) to be synthesized may be considered 
the same as a cascade compensator in the system. Thus by 
representing the problem in this manner the CALICO program 
may be used to approximate the transfer function cf the 
actual system by varying the parameters of this pseudo 

compensator (actually the transfer function of the system) 
in order tc minimize the cost function. The cost function, 
here represents the error between the measured frequency 
response data (entered as the desired profile) and the 
freguency response simulation calculated by the pregram. 
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Thus in effect the program minimizes the errcr between the 
transfer function approximating the system and the measured 
freguency response of the system at the discrete frequency 
points which form the profile input into the program. 

As is the case when using the program to design 
ccmpensa tcrs the user must assume a specific order fcr the 
numerator and denominator representing, in this case, the 
transfer function to be synthesized. While unfortunately 
the program will not vary the order of the numerator and 
denominator polynomials in order to obtain the best match to 
the measured freguency response, the user can employ the 
measured data in order to make an initial guess, as it were, 
of the form cf the transfer function. After this, one may 
scrutini2€ the results from the program, which graphs the 
measured and simulated responses, to make a decision if the 
form of the transfer function should be adjusted in order to 
more closely represent the system. Here also, a time domain 
simulation might be extremely useful in guiding decisions 
regarding the particular form cr modification that should be 
tried in crder to obtain a "good" model for the system under 
investigation. The fact that the user selects the specific 
orders cf the transfer function numerator and denominator 
and the minimization routine varies the parameters to cbtain 
a minimum cost function while leaving the orders unchanged 
presents the potential of obtaining a low crder model of a 
higher crder system. While for linear systems, with known 
transfer functions, there are numerous analytical techniques 
available for accomplishing this, the capability exists here 
of essentially obtaining a reduced order transfer function 
without a priori knowledge of the higher crder transfer 
function describing the system. The only data required by 
the program is the frequency response of the higher order 
system and nothing concerning the parameters of the high 
crder transfer function is required. 



128 



Since sene transfer functions may contain roots of the 
numerator polynomial in the right half of the s-plane 
(nenminiaum phase systems) this type of transfer function 
should cot te overlooked when attempting tc model a system. 
To this end a nonminimum phase option may be chosen ir the 
CALICO computer algorithm. When this option is selected the 
search area for the numerator coefficients is expanded to 
include both positive and negative values and the implicit 
constraints requiring left half plane roots for the 
numerator polynomial are eliminated. Selection of this 
option dees not guarantee that a nonminimum phase transfer 
function will result from the CALICO program but it does 
allow this type of transfer function to be considered as a 
feasible solution of the problem. 



B. TBANSI E B fONCIION SYNTHESIS: EXAMPLE PBCELEMS 



In this section several example problems are presented 
to illustrate the use of the CALICO program in determining 
transfer functions from frequency response data. These 
examples illustrate the validity of the basic concept of 
using the miniaizatipn algorithm to determine the transfer 
function, hut they are by no means an exhaustive 
presentation of all possible variations of cost functions 
and freguency ranges that may be employed in determining a 
transfer function representation of a system. For example, 
as will te shewn in the example problems, different transfer 
function representations may be obtained by considering 
different ranges of the measured frequency response. That 
is, while the order pf the numerator and denominator are 
fixed by the user, the parameter values returned from the 
program may vary depending upon the frequency range being 
considered. This nay be desirable if seperate low frequency 
and high frequency representations are being determined. 
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Also in using the type one cost function, which includes 
both phase and magnitude, the cost surface is different than 
if the magnitude difference alone were used in the cost 
function. This has an effect on the resulting parameter 
values, particularly when lower order transfer functions are 
chosen tc represent higher order systems. 

The necessary data and the input format employed are 
almost identical to those used when the program functicns in 
its originally designed compensator design mode. Figure V-1 
illustrates the input data deck used for the first synthesis 
example problem. As can be seen the only significant change 
is that the plant has been represented by a unity gain 
transfer function,. This has the effect that the resulting 
' , ccmpensatcr ,, represents the transfer function modeling the 
entire system. 
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Transfer Function Synthesis, Data Deck 
Figure V-l 



1 . 



Synthesis Example 

The freguency response data, shewn graphically in 
figure V-2, was generated from a known transfer function 
using a computer frequency response program. As can be 
seen, twerty magnitude and phase values were selected at 
discrete frequencies over the range from 0.1 radians/second 
to 30 radians/second. Examination of the freguency response 
data suggests a transfer function of the form of a corstant 
gain over a polynomial of some unknown order. This form 
will yield a constant magnitude at low frequencies and a 
decreasing magnitude as the frequency increases. Closer 
examination cf the magnitude curve in the range between 2 
and 20 radians/second indicates a slope of approximately 60 
cb/cecade, implying a polynomial of at least third order in 
the denominator. The phase curve extends from 10 degrees to 
330 degrees over the range cf frequencies considered. This 
indicates that the denominator may be of fourth order or 
higher. As a first guess a third order denominator is 
assumed. The results from the CALICO program are shown in 
figure ¥-2A through V-2H. The initial guess for the 
denominator iccts and the system gain were chosen as 1, as 
can be seen in figure V-2F. Figures V-2A and V-2D show that 
the resulting magnitude curve approximates the measured 
response well at the lower frequencies with a slight 
deviation begining to appear at the higher frequencies. The 
phase plot cf figure V-2B and the plot cf the difference 
between the actual phase and desired phase shewn in figure 
V-2E indicate a significant ameunt of deviation from the 
measured values in the higher frequency region. This 
suggests that another pole in the transfer function is 
necessary to mere accurately model the system. Thus, the 
program was again executed, this time using a fourth order 
deccminatcr. The results illustrated in figures V^2I 
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through V-2E indicate that this second form of the transfer 
furcticn acre accurately represents the system. Since the 
original transfer function is known in this case, an 
opcrtunity eiists to evaluate the performance of the program 
in terms of the accuracy of the transfer function returned 
after miniiii 2 ation of the cost function based on the input 
frequency response data. This of course would generally not 
be the case when trying to determine the transfer fucction 
of an actual system that was to be modeled. The parameter 
values returned from the CALICO program shown in figure V-20 
may be compared with the original transfer function given 
below. 

126.0 

T (s) = 

(s + 1) (s ♦ 2) <s + 5) (s + 10) 

As can be seen, t_he values of the transfer function roots 
returned from the program for the fourth order denominator 
are in agreement with the values of the original transfer 
furction used to generate the frequency response data. If 
for the system under consideration only the maguitude 
response were of interest, the transfer function consisting 
of a constant gain over a third order denominator might be 
of sufficient accuracy to model the system. This judgement, 
of course, rests with the user who is familiar with the 
intended use of the transfer function modeling the system. 
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20db 
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Magnitude vs. Frequency, Run #1 
Figure V-2A 



10.0 (rads/sec) 
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Figure V-2B 




Magnitude vs. Phase, Run #1 
Figure V-2C 
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Magnitude Difference vs. Frequency, Run #1 
Figure V-2D 
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Figure V-2E 
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TC r A L NU^HL-H CF TRIALS CALLliQ FUR = 10000 
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Computer Numerical Output, Run #1 
Figure V-2G 
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Computer Numerical Output, Run #1 
Figure V-2H 
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Magnitude vs. Frequency, Run #2 
Figure V-2I 
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Phase vs. Frequency 




Figure V-2K 
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Magnitude Difference vs. Frequency, Run #2 
Figure V-2L 
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Phase Difference vs. Frequency, Run #2 
Figure V-2M 
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Computer Numerical Output, Run #2 
Figure V-2P 



2 . 



Synthesis Exa mple 2 * 

As a second synthesis example, consider the 
frequency response data shown in figure V-3. Here the 
magnitude and phase response for a system are shown at 
twenty distinct freguency values over the range from 0.1 to 
100C.0 radians/seccnd. As can be seen from this figure the 
positive slope of 20 db/decade of the mgnitude curve at the 
lower frequencies suggests a transfer function with at least 
one zero. This is also confirmed by the phase response, 
which at the lower frequencies initially has a value of 90 
degrees. At the higher frequencies inspection of the 
magnitude plot shows a negative slope of approximately 40 
dt/decade, indicating the dominance of a third order 
denominator in this area. The phase plot approaching -180 
degrees at these higher frequencies also suggests this. 
Thus, initially a form for the transfer function consisting 
of a first order numerator and a third order denominator was 
selected to model the system from which the freguency 
response data was obtained. The results are shown in 
figures V-3A through V-3H. As can be seen from the 
graphical output presented in the figures, this particular 
form of transfer function appears to match the measured 
freguency response guite well over the entire range of 
frequencies considered. The roots and gain of the transfer 
function returned from the program are shown in figures V-3G 
and V-3H. As in the previous example it happens that the 
exact transfer function used to generate the measured data 
input to the CAIICC program is known. This transfer 
function is given below in order that the reader may compare 
the values of the original transfer function parameters with 
these returned from the program. 
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T(s) 



9000s 



(s ♦ 2.4) (s ♦ 5) (s ♦ 150) 

As can be seen from the data presented in figures V-3G and 
V-2H, the parameter values returned from the program are in 
close agreement with the original transfer function values. 
While it is unlikely that they will ever be exactly the same 
as the original values there is certainly net a significant 
difference between the two sets of values fer most 
engineering purposes. 
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Magnitude vs. Frequency 
Figure V-3A 



actual response 
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phase vs . Frequency 
Figure V-3B 




Figure V-3C 
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04db 
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Magnitude Difference vs. Frequency 
Figure V-3D 
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Phase Difference vs. Frequency 
Figure V-3E 
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Computer Numerical Output 
Figure V-3F 
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Computer Numerical Output 
Figure V-3G 



OPTIMIZED CCNPENSATCR TRANSFER FUNCTION DENCM I NATCR 
COEFFICIENTS IN ASCENDING POWERS CF S 
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Computer Numerical Output 
Figure V-3H 



3. 



Synthesis Example 3_. 

A freguency response of a system with resonant 
peaks, shcwr in figure V-4, is considered in this example. 
The same example is considered in ref. 19, where a 
ccmhinaticc cf graphical and analytical techniques are used 
to determine the transfer function. Initially a transfer 
function consisting of a first order numerator and second 
order denominator was assumed in trying to model the system. 
This choice was based on the low and high frequency 
magnitude curve slopes of approximately +20 db/decade and 
-20 db/decade respectively and the low and high frequency 
phase values cf approximately +90 degrees and -90 decrees. 
This form failed to match the dip in the magnitude curve in 
the region cf 10 radians/second and the jump in the phase 
curve in the same frequency area. After increasing the 
order of the numerator and denominator, a transfer furction 
consisting cf a third order numerator and fourth order 
denominator was found to provide a phase and gain curve that 
closely approximates the measured frequency response. The 
results cf this are shown in figures V-4A through V-4H. The 
numerical values of the transfer function parameters 
returned from the program closely approximate those of the 
actual system from which the frequency response measurements 
were obtained. The interested reader may verify this fact 
by comparing the values shown in figures V-4G and V-4E with 
these of the actual system given in ref. 19. 



162 




163 



Input Frequency Response Data 
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Magnitude vs. Frequency 
Figure V-4A 



actual response 
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Phase vs . Frequency 
Figure V-4B 




g -20db — 

Magnitude vs . Phase 
Figure V-4C 
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Magnitude Difference vs. Frequency 
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Computer Numerical Output 
Figure V-4H 



4 



Synthesis Example 4_. 

In this fiaial example problem a low order 

approximation to a known higher order system is illustrated. 
The original system transfer function under consideration 
consists of a seventh order numerator and an eight order 
denominator given by 

7 6 5 4 3 2 

16s + 463s +60 1 0s +36380s +122664s +222088s +185760s+40320 

8 7 6 5 4 3 2 

s +36s +546s +45 36s +22449s +67284s +118124s + 109584s + 40320 

An analytical methpd of approximating this high order 
transfer function is discussed in ref. 21. Initially 40 
discrete points, representing the frequency response cf the 
high order transfer function were selected over the 
frequency range from 0.01 radians/second tc 100.0 
radians/second. In prder to be able tc form a comparison 
with the reduced order model obtained in ref. 21, a transfer 
function consisting of a first order numerator and a second 
order deccirinator was assumed. The results returned from 
the program using this form of transfer function are shown 
in figures V-5A through V-5H. It is not surprising that 
there is a significant amount of deviation in the magnitude 
curve in the vicinity of the resonant effects of the higher 
order system nor tba<t the phase response cf this reduced 
order mcdel dees not match the higher order phase response 
in the higher frequency range. In order tc compare these 
results with analytical methods for forming reduced order 
models, the frequency response for the reduced order model 
of ref. 21 is shown in figures V-5I through V-5M. The 
numerical values returned from the CALICO pregram are cf the 
same order cf magnitude as those found by the analytical 
technique, but as can be seen frem the magnitude curves of 
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the two approaches (figuresV-5A and V-5I) , the analytical 
method matches the higher order transfer function extremely 
well at the lower and higher frequencies and makes little 
attempt tc match the higher order system in the ranee of 
frequencies at which the resenant peaks occur. The results 
returned frem the prpgram, on the other hand because they 
are a function cf the deviation of the response from the 
desired values over the entire range of frequencies, tend to 
distribute any error between the responses over the entire 
range of frequencies being considered. A slightly different 
representation of the system may be obtained by varying the 
freguency range under consideration and the type of cost 
function used in the program. For example, figures V-5N 
through V-50 illustrate the effects of using a type two cost 
function and limiting the freguency range from 0.1 to 100.0 
radians/seccnd. Again the decision as to which mcdel is 
best must be made by the user in light of the intended use 
of the recuced order representation. 
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Phase vs. Frequency, Run #1 
Figure V-5B 




Figure V-5C 
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Magnitude Difference vs. Frequency 
Figure V-5D 
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Figure V-5E 
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Computer Numerical Output, Run #1 
Figure V-5F 
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Computer Numerical Output 
Figure V-5H 



lOdb 
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Magnitude vs. Frequency, Analytical Transfer Function 
Figure V-5I 
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Phase vs. Frequency, Analytical Transfer Function 
Figure V-5J 




Figure V-5K 
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Magnitude Difference vs. Frequency, Analytical Transfer Function 

Figure V-5L 
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Phase Difference vs. Frequency, Analytical Transfer Function 

Figure V-5M 
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Magnitude vs. Frequency, Run #2 
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Phase vs. Frequency, Run #2 
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Magnitude vs. Phase, Run #2 
Figure V-5P 
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Magnitude Difference vs. Frequency f Run 
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Phase Difference vs. Frequency, Run #2 
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Computer Numerical Output, Run #2 
Figure V-5S 
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Computer Numerical Output, Run #2 
Figure V-5T 
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VI 



CONCLUSION HD SUGGESTED FUTURE ST UDIES 



The investigat ions carried out in this thesis illustrate 
the feasibility cf automating the classical design of 
compensators by specifying a set of desired values for the 
magnitude and phase response of the open loop system at a 
ouster cf discrete frequency values over the range of 
interest. The results achieved show that the complex 
airimizaticn method of M. J. Box is quite capable of 
handling the constrained minimization problem and 
conveniently provides a means of avoiding the numerical 
differentiation cf -the cost function, required with gradient 
search techniques. The incorporation cf stability criterion 
in the implicit constraints of the minimization algorithm 
also insures that the program will not design an unstable 
compensator. Ey specifying the desired response in terms of 
the magnitude and phase profiles, common frequency domain 
specif icaticns such as open loop bandwidth, phase margin, 
and gain margin have been incorporated into one common 
format and the necessity for specialized routines to check 
these specifications have been avoided. Those example 
problems presented for which known results had been 
established by other techniques show that the algorithm 
presented provides accuracy that is certainly within the 
necessary tolerances for most engineering applications of 
compensator design. 

Also, while not specifically designed for the purpose, 
the basic premise cf the algorithm's minimizing the 
difference between some desired frequency response and the 
frequency response of the system as calculated by the 
computer program while the transfer function parameters are 
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varied has been shewn to be effective in both determining a 
transfer function from measured frequency response data and 
in forming lew order models of known higher order linear 
systems. In the formulation of lew order models the choice 
of cost function and the frequency range specified may be 
usee tc emphasize different aspects of the freguency 
response which the low order model is tc approximate. 

While this work has shown that the basic algorithm 
presented can be used to effectively design series 
compensators based on frequency response specif ications, 
considerable work remains to be dene in this area. 
Specif ically , an option to select that only real poles and 
zercs be considered as viable solutions to the compensator 
design algorithm would not be a difficult addition tc the 
program and should prove guite useful. Also, the 
incorporation of a provision for handling pure time delays 
would not require significant modification of the program 
ceding and would extend the capability of the algorithm tc 
the solution of problems involving transport delays. 
Investigation of the method presented should be further 
pursued in terms of multiloop systems and systems involving 
compensation in the feedback loop. It may also prove 
worthwhile tc investigate extending this method to nonlinear 
components using describing function techniques. 

In the areas of transfer function synthesis and reduced 
order modeling, a more detailed investigation of both the 
effects of different cost functions and considering 
different frequency ranges should be conducted. In 
addition, it appears feasible to weight the cost furction 
over a particular frequency range by varying the density of 
the discrete frequencies at which the desired magnitude and 
phase prefiles are specified. Time did not permit a 
detailed investigation of this idea and further 
investigation into this area may be applicable to the 
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program utilization in both the areas of compensator 
and transfer function synthesis. 



design 
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APPENDIX A 



DESCRIPTION OF SUBROUTINES 



In this section a brief description and background of 
various subroutines is presented to aid the interested 
reader in better understanding the internal logic of the 
program. Also these descriptions are intended to facilitate 
the task of anyone desiring to make mcdif icatiors or 
additions to the program in order to handle systems cf a 
configuration ether than the unity feedback form assumed 
throughout this work. 
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SUBROUTINE CONORS 



The purpose cf this subroutine sub-program is to 
normalize tie coefficient cf the highest order term to 
unity. Ibis is accomplished in a simple, straight forward 
manner by dividing all the coefficients of the polynomial by 
the coefficient cf the highest order term. Thus a 
polynomial given by 

n 



P (S) = ^ a. 



i=0 

becomes cn return from this subroutine 

n 



P(s) 



<a 

■Z-r -- 1 



i=0 



where a is the coefficient cf the s term. 



The normalizing 



n 

factor a is also returned to the calling pregram, 
n 



It should be noted that the coefficients of the 
normalized pclynocial are returned under the same variable 
name as the input coefficients. 



Definitions cf the input and output parameters are: 

Input Variables 

V A real one-dimensional array containing the 

coefficients of the polynomial whose highest order 
coefficient is to be normalized to unity. Upon 
return from the subroutine this array will contain 
the normalized coefficients of the polynomial. 
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N The dimension cf the array containing the 

coefficients tc be normalized. 

VNOfiM The normalizing factor, returned tc the 

calling pregram, by which all coefficients have been 
divided. 
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SUBROUTINE PHASE 



This subroutine sub-program is used to calculate the 
phase angle cf a polynomial at a particular frequency. In 
order tc circumvent the difficulty associated with standard 
arctangent routines, which return values limited tc the 
range - v 0 > t , the phase contributions froa the 

individual factors are calculated and then summed to form 
the total phase anqle at a particular frequency. That is, 
if the polynomial is factored and represented as the product 
of first crder terms of the form 

n 

P(S) = YJ < s + O 

i=0 



where 


in general z 

i 


may be a 


complex number of the 


form 


z = 
i 


a + jfc , then 
i i 


at 


a 


particular frequency s = j 


<*» the 


phase 


contr itution 


of 


the 


individual factors may 


be 



evaluated as 

(b + ) 

i 

-1 

4> = tan 

i 

a 

i 

and the tctal phase angle at the particular frequency is the 
sum of these individual phase contributions. The phase 
contribution cf each of these first order terms will be 
within the principal values of the arctangent routine cf the 
computer and the problems generally associated with 
determination cf the phase for higher than first order terms 
are avoided. 
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Definitions of the input and output parameters are: 

Input Variables 

Bfi A real one- dimensional array containing the 

real farts cf the roots of the polynomial. 

fil A real one-dimensional array containing the 

imaginary parts of the roots of the polynomial. 

N An integer value specifying the number of 

rccts cf the polynomial. 

CHEGA The frequency in radians/sec at which the 

value cf the phase is to be computed. 

Output Variables 

TAN GD The total phase contribution of the 

polynomial given in degrees. 
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SUBROUTINE PINVBT 



The purpose cf this subroutine suh-prcgram is to 
rearrange tie coefficients of a polynomial in inverse order 
frcm that fed into the subroutine. Consider for example a 
polynomial arranged in descending powers of s 

n n- 1 

E (s) = a s ♦ a s +... +as+a 
n n-1 i o 

Upon return frcm this subroutine the coefficients would be 
arranged in ascending powers cf s. That is the polynomial 
would be arranged as fellows 

n- 1 n 

E (s) - a +as+...+a s +as 
oi n-1 n 

While both representations are mathematically equivalent 
the sequence cf the coefficients must be taken into account 
when using various polynomial manipulation techniques within 
the computer. The method used in reversing the sequence of 
the input coefficients consists of temporarily storing the 
first coefficient and sequentially shifting the remaining 
coefficients. The coefficient in temporary storage is then 
placed in the last cpefficient location. This process is 
continued with the use of nested do loops and each tine the 
shifting operation occurs the coefficient in temporary 
storage is placed one location lower than the previous 
number retrieved from temporary storage. 

This subroutine requires no ether subroutines or 
functions to carray put the procedure. The size of the one 
dimensional array of coefficients that may be rearranged is 
determined by the dimension cf the array in the calling 
program. 

Definiticns of the input and output parameters are: 
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Input Variables 



X (I) A real one dimensional arraj containing the 

coefficients to be rearranged. Upon return from the 
subroutine this array contains the input coefficients 
in reverse sequence from that entered. 

N The dimension of the array containing the 

coefficients to be rearranged. 
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SUBBOUTINE PMPY 



The purpose of this subroutine sub-program is to 
multiply the coefficients of one polynomial X (s) times the 
coefficients cf another polynomial Y (s) to obtain the 
coefficients cf the resulting polynomial Z (s) . That is, if 
it is desired to multiply two polynomials of the form 

2 n -1 n 

P (s) = a ♦ a s ♦ a s + ... + a s + as 

1 o i 2 n - 1 n 

and 

2 m -1 m 

P <s) = b + b s ♦ b s + ... ♦ b s + bs 

2 0i2 m- 1 m 

together in order to obtain a third resultant polynomial of 
the general fcrm 

2 i -1 i 

p (s) = c +CS+CS + ... + c s + c s 

3' 0 i 2 i-1 1 



this subroutine is used in accomplishing this. The order 
of the resulting polynomial is calculated by forming the sum 
of the orders cf the two input polynomials as 

i = n + m 

The various coefficients of the resulting polynomial, which 
are returned iE ascending powers of the independent 
variable, are calculated as the sums of products of the 
input coefficients by means of two nested do-loops. The 
arrangement cf the two do-loops is such that the appropriate 
product terms are summed in forming the respective 
coefficierts cf the resultant polynomial. All input and 
output coefficients must be real numbers. The size cf the 
polynomials that nay be multiplied are controlled by their 
respective dimension statements within the calling program. 



Definitions of the input and output parameters are: 
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Input Variables 



X A real one dimensional array containing the 

coefficients of the first polynomial tc be 
multiplied. These must be arranged in ascending 
powers of the independent variable. 

IDIMX The dimension of the array X c ontaining the 

coefficients pf the first polynomial. 

Y — — A real one dimensional array containing the 

coefficients of the second polynomial tc be 
multiplied. These must be arranged in ascending 
powers of the independent variable. 

IDIflY The dimension of the array I containing the 

coefficients of the second polynomial. 

Output Variables 

Z A real one dimensional array containing the 

coefficients of the resultant polynomial. These 
coefficients are arranged in ascending powers of the 
independent variable. 

IDIfiZ The calculated dimension of the array Z 

containing the coefficients of the resultant 
polynomial. 
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SUBBOOTINE PVAI 



The purpose cf this subroutine sub-prcgram is to 

evaluate a pclyncuial having the complex variable s = <r + j<*> 

and real coefficients a . For a given polynomial of order n 

i 

of the form 

n 

P (s) = ^ ^ a i* S 
i=0 

s is declared as a complex variable within the subroutine 
and the indicated operations of summation and multiplication 
are carried cut by a simple loop operation taking advantage 
cf the asscciative and distributive laws cf multiplication 
and addition. This subroutine requires nc other 
sub-programs tc perform its function. 

Definitions of the input and output parameters are: 

Input Variables 

A (I) A real one-dimensional array containing the 

coefficients cf the polynomial that is tc be 
evaluated. The coefficients are assumed to be 
arranged in ascending order according to the powers 
of s with zero explicitly input for the coefficient 
of any missing powers of s occuring in the 

polynomial. The maximum size of this array is 
controlled by the dimension of the array withir the 
main program. The values of the coefficients in this 
array remain unchanged after execution cf the 

subroutine. 
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NN The order of the polynomial tc he evaluated. 

This value is one less than the dimension of A (I) . 

PB The value of the real part cf the complex 

variable s at which the polynomial is to be 
evaluated. 

FI The value of the imaginary part cf the 

complex variable s at which the polynomial is to be 
evaluated. 

Output Variables 

V B — The real part of the complex resultant value 

of the polynomial evaluated at s = <r + j«>. 

VI The imaginary part of the complex resultant 

value cf the polynomial evaluated at s = <r + j*. 



208 



SUBROUTINE ROUTH 



The purpose cf this subroutine sub-program is to 
determine if a given input polynomial possesses right half 
plane roots by using the Routh stability criterion. The 
polynomial coefficients are input in descending powers of 
the independent variable. From these coefficients the Routh 
array is calculated and the first column of this array is 
searched for an indication of the existence of right half 
plane roots. Thus from an input polynomial cf the form 

n n- 1 

£ (s) = b s + b s + ... + bs + b 
n n-1 10 

the Routh array 
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n-4 


nd 

2 


row 


b 


b 


b 






n- 1 


n-3 


n-5 


rd 










3 


row 


a 

3 1 


a 

32 


a 

33 


th 










4 


row 


a 

4 1 


a 

42 


a 

43 



th 

n row 



st 

(n+1) row 



n 1 
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(n+1) 1 (n + 1) 2 (n+ 1 ) 3 
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is formed, where in general the entries from the third row 
onward are given by 

_ [ a {n-1> l]f (n-2) (i+1)] " f (n-2) l][ 3 (n-1) (i+1)] 

3 " ' ' ' ' " " — — .. . r m inrr nn — .. I- 

ci 

a (D-1) 1 

The first column of this array is then tested to determine 
if any of the elements are less than zero. Note that since 
all coefficients are limited to positive values by the 
search area restrictions in ether parts of the program no 
prevision has been made to handle the case cf the 

coefficients cf the polynomial being all negative. Thus a 
polynomial with all negative coefficients should be 

multiplied by -1.0 prior to using this particular 

subroutine fer testing for right half plane roots. 

Definitions of input and output parameters are: 

Input Variables 

Y A real one dimensional array containing the 

coefficients of the polynomial in descending powers 
of the independent variable. 

N The dimension of the array Y. 

Output Variables 

ISTABL Integer output indicating the absence or 

presence of right half plane roots 

Q--*-Indicates the polynomial has no right half 
plane rpots. 

1 Indicates the polynomial has right half 

plane roots. 
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SUBROUTINE SEMEL 



The purpose of this subroutine sub-program is to 
calculate the coefficients of a polynomial from the rccts of 
the polynomial. The calculation is done using complex 
arithmetic in order to be able to compute the polynomial 
coefficients when complex roots are present. The resulting 
polynomial coefficients are, however assumed to be real. 
Thus, if complex roots are present they must be in conjugate 
pairs. ilsc the coefficient of the highest order term is 
set egual to unity. Consider a polynomial of degree n 
having n roots. This may be expressed as 



where in general the z 's may be complex. The coefficients 

i 

of the polynomial that results when the above indicated 
multiplication is performed may be computed by the formula 



n 




(s + z.) 
1 



i=0 



D-l r — > 

A. = (-1) • / (product of roots taken n-i at a time) 



i 



Thus the coefficients of the polynomial are given by 
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z «z •z »z • . . 

12 3 4 



a 

o 



. «z 

n 



a = z «z «z •...•z + z •z «z •...•z 

1 12 3 n-1 2 3 4 n 



a = z +z +z ♦ z +z ♦ . . . + z 

n- 1 1 2 3 4 5 n 

a = 1 

Q 

The largest crdet polynomial that this particular subroutine 

th 

will handle is a 2-0 order. If larger polynomials are to 

be considered the dimension of R within the subroutine must 
be changed. 

Definitions of the input and output parameters are: 

Input Variables 



BE A real one-dimensional array containing the 

real parts of the polynomial roots. 

BI A real one- dimensional array containing the 

imaginary parts of the polynomial rccts. 

K An integer specifying the order of the 

polynomial or equivalently the dimension of the 
arrays containing the real and imaginary parts cf the 
rccts . 



Cf 



Output Variables 

A real one- dimensional array containing the 
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resulting coefficients cf the polynomial arranged 
ascending powers of the independent variable, 
dimension of this array must be one larger than 
dimensions of the arrays containing the real 
imaginary parts of the roots. 
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APPENDIX B 



PROGBAM LISTING 
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APPENDIX C 



DESCRIPTION OF FOBTRAN VARIAELE NAMES 
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FOBTBAN 



IABIABLI 


OMIT 


DESCRIPTION 


A (I) 


None 


Array containing the coefficients 
of the plant transfer fuEction 
numerator . 


ADAMG (I) 


db 


Array containing the desired cpen 
loop magnitude response in db at 
the discrete frequencies specified 
by the designer. 


AMAG (I) 


db 


Output array containing the 
magnitude of the open loop 
frequency response in db at the 
discrete frequencies specified. 


ANGDEE 


Degrees 


Phase angle contribution of the 
open loop system transfer function 
denominator at a specific 
frequency. 


ANGDEB 


Badians 


Phase angle contribution of the 
open loop system transfer function 
denominator at a specific 
frequency . 


ANGNUE 


Degrees 


Phase angle contribution of the 
open loop system transfer function 
numerator at a specific frequency. 


ANGNOB 


Badians 


Phase angle contribution of the 
open loop system transfer function 
numerator at a specific frequency. 


E (I) 


None 


Array containing the coefficients 
of the plant transfer function 
denominator . 


CO (I) 


None 


One dimensional working array used 
primarily in reading polynomial 
coefficients and inputing data to 
root finding subroutines. 


COF (I) 


None 


One dimensional working array used 
primarily to store resulting 
coefficients used in checking 
accuracy of root finding 
subroutines. 


COFASN (I) 


None 


Array containing the coefficients 
of the cpen loop transfer function 
numerator. 
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FORTRAN 

VARIABLE OMIT DESCRIPTION 



COFASD (I) 


None 


Array containing the coefficients 
of the cpen loop system transfer 
function denominator . 


COFWRK (I) 


None 


One dimensional working array used 
primarily to temporarily store 
polynomial coefficients in 
performing polynomial 
manipulations. 


DELPHI 


Degrees 


Difference between two numerator 
phase angle contributions at 
successive discrete freguencies. 


DELPSI 


Degrees 


Difference between two denominator 
phase angle contributions at 
successive discrete freguencies. 


DIFF1 (I) 


d b 


Difference between the desired and 
actual magnitude of the open loop 
freguency response at the discrete 
frequency values specified in the 
input. 


DIFF2 (I) 


Degrees 


Difference between the desired and 
actual phase of the open lccp 
freguency response at the discrete 
freguency values specified in the 
input. 


DM AG (I) 


None 


Array containing the desired 
magnitude response at the discrete 
freguencies specified by the 
designer. 


DQPHSE (I) 


Degrees 


Input array containing the desired 
phase profile of the open loop 
system at the discrete freguercies 
cnosen. 


EWFQ 


Radians 


Interval to be used if discrete 
freguency values are to be 
incremented linearly. 


FREQ (I) 


None 


Logarithm of the discrete 
freguencies used in constructing 
abscissa for Bode plots. 


GAIN 


None 


Gain of plant transfer function. 


GAINCC 


None 


Gain of compensator transfer 
function. 



295 



FCBTBAN 

¥ARIABIE 


ONIT 


DESCRIPTION 


ICB 


None 


Input integer to be set egual to 
unity if the desired open loop 
magnitude response is input in db. 


IEBBXR 


None 


Output integer indicating the 
completion code of the minimization 
routine. 


INFORM 


None 


Alphabetic input indicating if 
coefficients are tc be read in 
factored or polynomial form. 


IPIOT 


None 


Integer input, to be set to UEity 
if CAICOMP plots of results are 
desired. 


IZOH 


None 


Integer input set egual to unity if 
there is a zero order hold in tne 
error signal path prior to the 
plant ana compensator. 


KNOW 


None 


Integer input indicating if 
discrete frequency values are to be 
read in or incremented linearly 
from the minimum freguency. 


NEODE 


None 


Integer input to be set equal to 
unity if Bode plot is not desired. 


NCD 


None 


Order of compensator transfer 
function denominator. 


NCN 


None 


Order of compensator transfer 
function numerator. 


NICHCI 


None 


Integer input to be set equal tc 
unity if Nichol*s plot is not 
desired. 


NMINFS 


None 


Integer input to be set equal to 
unity if a nonminimum phase 
solution is to be allowed. 


NOMEG 


None 


Integer specifying the total number 
of discrete frequency points being 
considered in the input freguency 
response profile. This number must 
be { 500. 


NSD 


None 


Order fo the system open lccp 
transfer function denominator. 



FCBTBAN 



ABI ABLE 


UNIT 


DESCRIPTION 


MSN 


None 


Order of the system open loop 
transfer function numerator. 


NTFD 


None 


Order of plant transfer function 
denominator. 


NTFN 


None 


Order of plant transfer function 
numerator. 


NIQST 


None 


Integer input to be set equal tc 
unity if Nyguist plot is not 
desired. 


M 3 


None 


Dimension of array containing plant 
transfer function numerator 
coefficients. 


N4 


None 


Dimension of array containginc 
plant transfer function denominator 
coefficients. 


CDENI 


None 


Value of the immaginary part of the 
system open loop transfer function 
denominator evaluated at a 
particular frequency. 


CDENR 


None 


Value of the real part of the 
system open loop transfer function 
denominator evaluated at a 
particular frequency. 


CM AG (I) 


None 


Output array containing the 
magnitude of the frequency response 
computed for the open loop system 
transfer function at the discrete 
frequencies specified. 


CNOMI 


None 


Value of the immaqinary part cf the 
system open loop transfer function 
numerator evaluated at a particular 
frequency. 


CNUMR 


None 


Value of the real part of the 
system open loop transfer function 
numerator evaluated at a particular 
frequency . 


CPHSD (I) 


Degrees 


Output array containing the phase 
of the frequency response of the 
open loop system at the discrete 
frequencies specified. 
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FORTRAN 

VARIABLE 



OMIT 



DESCRIPTION 



CPHSR (I) Radians Output array containing the phase 

of the frequency response of the 
open loop system transfer function 
at the discrete frequencies 
specified. 



SAMAG (I) 


db 


Array containing the magnitude of 
the simulation of the resulting 
optimized system transfer function 


SAPHSE (I) 


Degrees 


Array containing the phase cf the 
simulation of the resulting 
optimized system transfer function 


1 


Secs 


Sampling period if a zero order 
hold is present in the system. 


WFREq<I) 


Radians 


Array containing the discrete 
frequencies specified by the 
designer. 


KMAX 


Radians 


Input specifying the maximum 
frequency value of the range cf 
frequencies being considered. 


REIN 


Radians 


Input specifying the minimum 
frequency value of the range cf 
frequencies being considered. 


NS 


Radians 


Output giving the sampling 
frequency if a zero order hold is 
present m the system. 


XS(I) 


None 


Array in which values to be varied 



by the minimizaticn routine are 
stored. At the start this array 
contains the initial guess of the 
coefficient values and at 
completion the array contains the 
coefficient values that yield the 
minimum cost function. 



298 



LIST OF REFERENCES 



Massachusetts Institute of Technology, Electronic 
Systems laboratory Report ESL-R-409, On the Design of a 
Digital Ccmpu ter Program for t he Design cT'Fee'dEack 
C*3m censatcrs In Transfer Function form# by "07 Itlians, 
T3ec€lE€l*77597 



Eeveridce, G. S. G. and Schechter, R. S.* 

Optimization: Theory and Practice# p.453-456, 

T3cT3fam*RITI7 197TJ7 — 

Bode, H. W., Anal ysis and Fee dbac k Amplif ier Design, D. 
Van Kostrand Company, 7345. 

Box, M. 0., 11 A ;New Method of Constrained Optimization 
and a Comparison With Other Methods", Computer Journal, 
p. 42-52 , 8 April 1965. 

Aerospace Corporation Report TR- 1001 (2307) -20 , The 
Application of Modern Computing Technology to Control 
System Inalysis and Des ign ^rdPl ems # ~Ey~T7 C7 CcfleyT 
3une*79E77 * 



Jaccfcy, S. L. S. , Kovalik# J. S., and Pizzo, J. T., 
Iterative Methods for Nonlinear Optimization Problems, 
1773*5=737 ,**Frentice=HaII7~7372. 



Kuo * F. f., and Kaiser, J. F., Editors. S yste ms 
Analy sis by Digital Compu ter , Wiley, 1 96o . * 

Leondes, C. T., Editor, Control and Dynamic Systems 
Advances in Theory and Applic atio ns# v. 77, p. 2*5-743, 
Icadellc Press, 73747 ** 

Melsa# J. L. and Jones, S. K.. Computer Programs for 
Computational Assistance in the Study cl l ine ar Control 
Theory, 237 ed . , McCfaw Hill, *7973. * * * 

Melsa, J. L. and Schultz, D. G., Linear C ont rol 
Systems, p. 513-591, McGraw Hill, *73337 



NASA Contractor Report CR-2248, An Innovative Approach 
to Compensator Design, by Jerrell**E. HitcleTI and 
*Billie I7**McDanieI7 Jr. , May 1973. 



Oldentcrger, R- , Editor, Fre que ncy Respo nse # The 
MacMillian Co., 1956. 



Page, J. A., Automated Design of Feedback Control 
Systems, Fh.D .“Thesis , University of*CaIif ornia , los 
Xng eles , 19 70. 



Rosenbrock, B. H. , Computer Aided Control System 
Design# Academic Press, 7374. * 

Sanathanan, C. K. and Tsukui, H., "Synthesis of 
Transfer Function from Frequency Response Data", 
Ihternaticnal Journal of Systems Science# v. 5, no. 4, 

p747=*547*3ay“ 13747 



299 



IISI OF REFERENCES (cont.) 



16 . 

17 . 

18 . 

19 . 

20 . 

21 . 




Thaler, G. J., D esig p of F eedbac k Sy st ems, Dowden, 
Hutchinscn, and "Boss, 1373. “ 



Thaler. G. J. and Brown, B. G. . Analysis and Design of 
Feedb ack Control Syste ms, 2d. ed. , HcGraw Hill, TS60 . 



Truxal, J. G., Introductory System Enqineerinq, FcGraw 
Hill, 1972. 



Office cf Naval Besearch London Conference Report 
C-20-73, Conference on "Computer Aided Control System 
Desig n" , 'Ey 1. H7~UaynicTc, 10 SepTem'Ber“l'973. 



Computer Aided Control System Design Conference 
Publication Number 96. A Method for Linear Systems Order 
Seduction, by S- C. Cnuang, p7213-2T8,~ 1 Ipril 1973. 



300 



INITIAL DISTRIBUTION LIST 



Nc. Copies 

1. Defense Eocu mentation Center 2 

Caaercc Station 

Alexandria, Virginia 22314 

2. Library, Code 0212 2 

Naval Ecstgraduate School 

Monterey, California 93940 

3. Department Chairman, Code 62 2 

Naval Ecstgraduate School 

Monterey, California 93940 

4. Professor Gecrge J. Thaler, Code 62TR 5 

Naval Ecstgraduate School 

Monterey, California 93940 

5. Lt. Anthony J. Mancini, USN 3 

1514 Electric St. 

Scrantcn, Pa. J8509 

6. Dr. Jerrel R. Mitchell 1 

Department of Electrical Engineering 

Mississippi State University 
Mississippi State, Mississippi 39762 

7. Mr. Jay Cameron 1 

IBM 

Monterey and Cottle Roads 
San Jcse , California 95114 



301 



8 



1 



Mr. Peter Gramata 
IBM 

Monterey and Cottle Roads 
San Jose, California 95114 

9. Mr. Ectert McIntosh 1 

IBM 

Monterey and Cottle Roads 
San Jose, California 95114 

10. Mr. Martin Dost 1 

IBM 

Monterey and Cottle Roads 
San Jose, California 95114 

11. Mr. Ben Palmer 1 

IBM 

Monterey and Cottle Roads 
San Jose, California 95114 

12. Mr. S. M . Syn 1 

IBM 

Monterey and Cottle Roads 
San Jose, California 95114 

13. Professor A. G. J. MacFarlane 1 

Dept, cf Electrical Engineering 

University of Manchester 
Inst, cf science and Technology 
Manchester England 067609 

14. Dr. Eui-Tien Rung 1 

Dept, cf Applied Sciences 

Universite du Quebec a Chicoutimi 
930 est, rue Jacgues-Cartier 
Chicoutimi, Quebec 
G7H 2E1 



302 



15 



1 



Mr. Dave Noble 
IBM 

Monterey and Cottle Roads 
San Jose, California 95114 



303 



Thesis 
M2 7 86 
c.l 



166101 



Mancini . 

Computer aided. 

system desi9 n ^ s ' 
frequency domain 
f ications. 



control 

ng 

speci - 



Thesis 
M2 786 
c.l 



. . 166101 

Mane i ni 

Computer aided contro’ 
system design using ( 

frequency domain speci- 1 
f ications. 



